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Abstract 


This paper studies the relations among system parameters, uniqueness, and stability 
of equilibria, for kinetic systems given in the form of polynomial ODEs. Such models are 
commonly used to describe the dynamics of nonnegative systems, with a wide range of 
application fields such as chemistry, systems biology, process modeling or even transportation 
systems. Using a flux-based description of kinetic models, a canonical representation of the 
set of all possible feasible equilibria is developed. 

The characterization is made in terms of strictly stable compartmental matrices to define 
the so-called family of solutions. Feasibility is imposed by a set of constraints, which are linear 
on a log-transformed space of complexes, and relate to the kernel of a matrix, the columns 
of which span the stoichiometric subspace. One particularly interesting representation of 
these constraints can be expressed in terms of a class of monotonous decreasing functions. 
This allows connections to be established with classical results in CRNT that relate to the 
existence and uniqueness of equilibria along positive stoichiometric compatibility classes. 

In particular, monotonicity can be employed to identify regions in the set of possible 
reaction rate coefficients leading to complex balancing, and to conclude uniqueness of equi¬ 
libria for a class of positive deficiency networks. The latter result might support constructing 
an alternative proof of the well-known deficiency one theorem. The developed notions and 
results are illustrated through examples. 
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Description 


n-dimensional real space 

n-dimensional positive (resp. negative) orthant 

n-dimensional non-negative (resp. non-positive) orthant 

each element of the vector x is positive (resp. negative) 

each element of the vector x is non-negative (resp. non-positive) 

The n-dimensional vector with each element being one 

the ith standard basis vector in R" 

diagonal matrix X>(x) £ R"’^" with components of x in the diagonal 
number of species 
number of complexes 

number irreversible chemical reaction steps in the network 

rate of the reaction from complex i to complex j 

rate coefficient of the reaction from complex i to complex j 

the number of linkage classes 

integer for indexing linkage classes 

the Ath linkage class 

the number of complexes in linkage class no. A 

index of the reference complex in linkage class no. A 

vector of concentrations (state variables) 

m X n dimensional molecularity matrix 

monomial function of the kinetic dynamics, 4>i{c) = YljLi 

net reaction rate corresponding to complex i 

subspace containing im{Ak) 

stoichiometric subspace 

dimension of the stoichiometric subspace 

deficiency of the reaction network 


Defining/introducing eqn. 
or (sub)section 


sec. 2, sec.3 
eq. ([28ll.(|73ll 
sec. 2 
sec. 2 
sec. 2 
eq. dU 
eq. dU 
sec. 2 
sec. 2 
sec. 2 

sec. 2 (footnote 1) 

subsec. 2.1 

sec. 2 

eq. dH) 

eq. dB) 

eq. dS]) 

eq. (fTOl) 

eq. (fTTI) 

subsec. 2.2 

eq. (fT^ 


1 Introduction 


Deterministic reaction networks obeying mass action law (MAL) kinetics form an important 
subclass of kinetic systems, which in spite of their apparent simplicity, are able to describe a 
rich variety of dynamical behavior, that includes multiple equilibria conditions, oscillations or 
even chaos [isiie]. Such networks are typically employed to describe the dynamics of open or 
closed chemical reaction systems, but over the last years they proved useful in modeling other 
system classes as well. 

Reaction networks belong to the class of nonnegative (or positive) systems, the main char¬ 
acteristics of which being that the non-negative orthant is invariant for the dynamics. The 
application field of nonnegative systems extends far beyond chemistry, and includes dynamical 
models whose state variables are naturally nonnegative, as it is the case of biological systems 
in their many scales (from cells to ecological systems), or systems that can be transformed 
to be nonnegative, such as certain process models (e.g. heat exchangers, distillation columns, 
convection networks), economic, transportation or stochastic models [20]. With an appropriate 
selection of coordinates, even many classical mechanical and electrical models can be described 
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in the nonnegative framework. 


The main specialty of reaction networks within nonnegative polynomial models is the lack of so- 
called cross-effects, what defines an additional constraint between the monomial coefficients and 
exponents m- Still, the class of reaction network models is quite wide, and many non-chemical 
models can be brought into a kinetic form using simple transformations PSH]. Widely used 
examples of kinetic systems are compartmental systems m and Lotka-Volterra models into 
which most smooth nonlinear ODEs can be embedded [33]. These facts, clearly underline the 
importance of reaction network models and motivate us to attempt to look at general dynamical 
models through the glasses of kinetic systems. 

The study of the relationships between chemical reaction structure and dynamic behaviour is 
the purpose of Chemical Reaction Network Theory (CRNT), a program formally proposed and 
developed in la 0 HOj. One of the earliest results on the relation between the solutions of 
nonlinear dynamical systems (including kinetic systems) and their associated directed graphs 
is published in [52] . Important cycle-related conditions on the stability of kinetic systems were 
given in (TOj. An extensive stability analysis of reaction networks using algebraic and graph- 
theoretical tools can be found in m- Thermodynamically motivated Lyapunov-function-based 
stability analysis of kinetic systems, considering certain frequently applied model-simplification 
steps, is proposed in [32] . 

The seminal works in [351121] (collected in their most comprehensive form in [23] 1 explored the 
dynamic properties of MAL complex chemical systems, and contributed to equip CRNT with 
a mathematical formalism that has prevailed to present. It is important to remark here, that 
several different network structures may correspond to the same kinetic differential equations m 
I18j . Therefore, important network properties such as deficiency, weak reversibility or complex 
balancing, may vary among the possible reaction structures belonging to the same ODE model 
(see, e.g. (SHIES])- 

One fundamental problem in CRNT is to decide from the structure or parameters of the network, 
whether it can exhibit or not multiple equilibria. An important early result in this field is the 
rigorous proof of the existence and uniqueness of thermodynamic equilibrium in a mixture of 
chemically reacting ideal gases [53]. The motivation in [49] was the computational analysis of 
large thermodynamical models. The work contains fundamental results about the existence and 
uniqueness of compositions minimizing the free energy. 

In answering the questions about the properties of equilibria, the concept of network deficiency 
(a number that relates to reaction network structure and stoichiometry) has become central to 
characterize the network behavior. Two essential results of CRNT are the well-known deficiency 
zero and deficiency one theorems |23l [23] which (besides other important results) establish 
conditions for networks to have exactly one equilibrium point in each positive stoichiometric 
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compatibility class |23j . This suggests network robustness with respect to parameter variability, 
and underlines the importance of the kinetic system class in general nonlinear systems theory. 

CRNT has received renewed interest over the last years, particularly in the area of systems 
biology, because of its potential to explore and to analyse complex behavior and functionality 
in biological systems (e.g. [laiiaiaiis]). Most efforts were dedicated to investigate the re¬ 
lationships between reaction network structure and dynamic behaviour. In this regard, special 
mention should be made of the so-called injectivity property, investigated as a condition that re¬ 
lates to the singularity (or not) of the determinant of the Jacobian associated to a given dynamic 
system [16]. Algebraic and graph theoretical methods have been devised to check injectivity, 
and therefore uniqueness of equilibria [13 Ej. In the same direction, extensions to cope with 
instabilities have been developed in [l2|- From different perspectives, a number of necessary and 
sufficient conditions for a given network structure and stoichiometry to accommodate multiple 
equilibria have been also recently proposed in [mss]. 

A particularly interesting class of chemical networks are the reversible ones, either in the strict 
thermodynamic sense, in which every elementary reaction step is reversible, or in a weak re¬ 
versibility context. Reversibility leads to a particular set of positive equilibria which is known 
as detailed balance if each reaction step is equilibrated by a reverse one, or complex balanced if 
the network is weakly reversible. 

At this point, it must be remarked that equilibrium should be understood along the sequel in 
the sense given in dynamic systems, irrespectively of whether it corresponds to thermodynamic 
equilibrium or to a particular steady-state on a chemical reactor. Note, however, that in agre- 
ment with thermodynamics, instabilities in the dynamics of reaction systems (when taking place 
on a homogeneous medium in isothermal conditions) require the reaction domain to be open to 
mass exchange with the environment. 

Because of microreversibility, most chemical systems, when closed to mass and energy exchanges 
with the environment, satisfy the principle of detailed balance equilibrium, resulting into stable 
equilibria [29] • As discussed in m and [28], irreversibility can be allowed within a reaction 
network, as limit cases of reversible steps under a thermodynamic consistency condition (known 
as the Wegscheider condition) which necessarily assumes microreversibility. 

The notion of complex balancing (also known as cyclic balancing or semi-detailed balancing), 
on the other hand, generalizes the detailed balance condition to any weakly reversible network. 
The structure of complex balanced systems has been explored in m and shown to be a toric 
variety with unique and stable equilibrium points (see also m)- Extensions to cope with more 
general classes of kinetic systems have been investigated in |43l I46j . 

It is important to mention here the recent fundamental results on the proof of the Global 
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Attractor Conjecture which says that any equilibrium point of a complex balanced mass action 
system is globally stable. A proof for the single linkage class case was given in [3], while a 
possible general proof based on differential inclusions was described in [Hj • 

CRNT, as it stands nowadays within the field of applied mathematics, offers an extraordinary 
potential in system’s theory for analysis and design of complex dynamic systems of polynomial 
type, what in turn may cover a wide spectrum of chemical and biological systems. Unfortunately, 
many of its results remain at a large extent unexploited, when not unnoticed, in the fields of 
process systems and engineering. 

Among the reasons that hamper application might be certain advanced mathematical tools and 
the intensive use of graph theory that are often not well-enough known to engineers. Some 
practical questions that demand attention relate to the link between dynamic behavior of a 
given mechanism and parameter sets (reaction rate coefficients), or to the design of a chemi¬ 
cal/biochemical network with some pre-specified behaviour (e.g bistable, oscillatory, etc). 

In this contribution we present some conditions that ensure feasibility of equilibrium solutions 
for weakly reversible mass action law (MAL) systems. They are linked to the notion of family of 
solutions, a concept originally derived in [HI SS] to study multiplicity phenomena as a function 
of network parameters. 

In deriving what it will be referred in the sequel as feasibility conditions, we exploit a flux- 
based form of the model equation. Within such structure, the time evolution of the species 
concentration vector is expressed as the product of a matrix denoted by S, whose columns span 
the stoichiometric subspace of the reaction system, and a vector function that is related to 
concentrations through a class of stable Metzler matrices [6]. 

As we will show, feasibility relates to the orthogonality between a log-transformed vector func¬ 
tion of reaction complexes and the kernel of matrix S. Based on this observation, feasibility 
conditions will be expressed in terms of certain functions that can be employed to identify admis¬ 
sible equilibria within the positive orthant of the concentration space. It will be shown that such 
functions are monotonous in their respective argument and take the zero within their domain, 
what will allow us to establish links with existence and uniqueness of equilibria along positive 
stoichiometric compatibility classes for MAL kinetic systems. In this context, connections be¬ 
tween monotonicity and two classical results in CRNT theory that relate to complex balanced 
equilibrium [MIES], and to a class of positive deficiency networks [22l[23l[21], will be discussed. 

Finally, it must be remarked that the potential interest of the notion of complex balancing in 
the context of process control is to characterize stable operation regimes in open systems, where 
the principle of detailed balance does not necessarily hold. This may allow, for instance, the 
selection or manipulation of exchange fluxes so to preserve stability of the resulting (open to the 
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environment) reaction system, via appropriate process optimization and/or feed-back control 
(see e.g. m)- Future directions may also involve the detection or design of networks having 
multiple equilibria. 

The paper is organized as follows: Section [2] introduces a formal description of chemical reaction 
networks. The graph structure underlying a reaction network, and its algebraic counterpart, 
will be described in Section [3l Section 0] presents a flux-based form canonical representation of 
the equilibrium set, that includes some feasibility conditions. Relationships between network 
structure and monotonicity of feasibility conditions will be established in Section [5j Connec¬ 
tions between monotonicity of feasibility functions and some classical results on uniqueness and 
stability of equilibria will be discussed in Section [UJ 


2 Preliminaries: Reaction Network Structure and Dynamics 


Let m be the number of chemical species which react by p irreversible chemical reaction steps, 
and c € the corresponding vector of species concentrations, defined as mole number per unit 
of volume. Each reaction step transforms some set of chemicals, usually referred to as reactants, 
into a set of reaction products. In CRNT, reactants and reaction products receive the name of 
reaction complexes. Complexes and reaction steps describe a graph where complexes correspond 
to nodes and reaction steps to directed edges. 

Formally, a graph involving n complexes {Ci,... ,Cn} linked by irreversible reaction steps can be 
constructed by associating to each complex i a set 2j with n integer elements, and a vector 
The elements of the set Xj are the indices of the complexes that are directly reachable (i.e. by 
one reaction step) from Cj. From now on, we will refer to each complex Ci by the corresponding 
index i. Vector y* G M”* has as entries the (positive) stoichiometric coefficients of the molecular 
species that participate in complex i. 

The graph structure is then built by linking every complex i to j G Xj. This process results in 
a number i of connected components known in CRNT as linkage classes. For each linkage class 
A = we define the set Cx which contains as elements the indexes of the complexes that 

belong to that linkage clasji]. 

Complexes are connected within a linkage class by sequences of irreversible reaction steps that 
define directed paths. Two complexes are strongly linked if they can be mutually reached from 
each other by directed paths (trivially, every complex is strongly linked to itself). A maximal set 
of pairwise strongly linked complexes dehnes a strong terminal linkage class if no other complex 

^To be precise, the set Cx is that containing as elements Cx = {*i, * 2 ,..., ijv;,,}, with Nx = JV{£,x), being ij the 
cardinality associated to complex Ci^, and Af{-) the operator which indicates the number of elements in the set. 


6 



can be reached from its nodes. In this work we will consider only networks in which every linkage 
class contains just one strong terminal linkage class. 

A linkage class C\ is said to be weakly reversible if any pair of its complexes is strongly linked. 
Weakly reversible networks are those composed by weakly reversible linkage classes. A particular 
type of weakly reversible linkage class is a reversible linkage class if each reaction step is itself 
reversible, so that for every i and j € Xj, we have that i G Xj. The rate Rij, at which a set of 
reactants in complex i is transformed into a set of products in complex j, will be assumed to be 
mass action, so that: 

m 

Rij{c) = kijtpiic), with V'j(c) = WcY = (1) 

i=i 

where yi is the stoichiometric vector corresponding to complex C,. The reaction systems we 
consider in this work will take place under isothermal conditions, what makes any reaction rate 
parameter kij{> 0) constant. Whenever c is a strictly positive vector, the following alternative 
representation for ipi{c) may be more convenient: 

\mpi{c) = yjlnc, (2) 

where the natural logarithm operator In(-) acts on any vector element-wise. Let ip : M^q —>■ M”q 
be the vector containing as entries the monomials described in ([1]), then the previous expression 
can be written in matrix form as: 

In 1 / 7 ( 0 ) = F"*" Inc, (3) 

where Y G jg so-called molecularity matrix which collects as columns the stoichiometric 

vectors y^ G associated to the complexes of the network. 

2.1 The dynamics of reaction networks 

Following the classical work by Feinberg |21] . the time evolution of species concentrations on 
a well-mixed reaction medium at constant temperature can be described by a set of ordinary 
differential equations that we write as: 

c = Y ■ Ak{'ip{c)) = Y (4) 

A 

where A^ : M” —M” is a linear operator dehned as: 

= X] (5) 

i&Cx jeii 

with Bi G M"' denoting the ith standard unit vector employed to represent axes on a cartesian 
coordinate system. Let us dehne the net reaction rate flux around a complex i, as the signed 
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sum of in- and out-flowing fluxes, i.e. as a function (j)i : 


^>0 


of the form: 


{j\ieXj} j£li 


(6) 


where the first summation at the right hand side extends to all source complexes j in the network 
from which there exists a reaction step to product complex i, and is represented by {j\i G Ij}. 

We can express in ([5]) in terms of fluxes ([6]), by selecting any reference complex j\ G 

and adding and subtracting Sj^ from the right hand side of ([5]) so that: 

iGCx j€li ieCx \j&Ii / 

After switching sub indexes, re-ordering the summations for the first term at the right hand side 
and making use of ([6]), we get the following equivalent expressions: 




AtW = X] - ^3 

ieCx \{i|iex,} / \j&Xi / 

= (7) 


i€Cx 


For convenience, the reference complex will be chosen from the corresponding strong terminal 
linkage class. Since vectors are orthogonal, by using ([7]), we have that = ejA^{rp) for 

every i £ Cx. Let ux = '^i^Cx have that ui^A^{il3) = 0 and therefore: 

<PiW = ( = 0 - ( 8 ) 


i^Cx 


<i€Cx 


Note that fluxes in ([6]) (as well as the linear operator in ([5|)) are implicitly dependent on the 
reaction rate coefficients associated to the reaction steps in the linkage class. By inspection of 
it can be concluded that the image of Afc('0) lies on the subspace A defined as follows: 


A = Ai + • • • + Aa + • • • + A^, 


where 


Aa = spanjej — \ i G Cx} for A = 1, • • 

and the sum of vector spaces Vi and V 2 is defined as: 

V1 + V2 = {ui + U 2 I Ui G w, V 2 ^ V2}. 


(9) 

( 10 ) 


Since vectors in {si — \ i £ Cx} are linearly independent, they form a basis for the subspace 

Aa, thus dim(AA) = Nx — 1. In addition, since the subspaces Aa are orthogonal: 

dim(A) = ^(A^a — 1) = n — 1. 

A 
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This implies that Ak{il)) = 0 if and only if = 0 for all i G Ua/1a- Consequently, if a 

positive concentration vector c exists compatible with a zero flux condition for every complex in 
the network, that vector should be an equilibrium for system ([3]). Such equilibrium condition, 
known as complex balanced |34] . is formally dehned as follows: 

Definition 2.1 (Complex Balanced Eqnilibrium) Any vector c* > 0 such that (/>i(t/>(c*)) = 
0 (Eqn for every i = 1, ■ ■ ■ ,n is called a complex balanced equilibrium solution. 

A subclass of complex balanced equilibrium, particularly meaningful from a thermodynamic 
point of view as it relates to microreversibility f |37l [29]), is the detailed balance equilibrium 
which we define next: 

Definition 2.2 (Detailed Balance Eqnilibrium) If the network is reversible (i.e. for every 
i and j £ li, we have that i £ Ij) any vector c* > 0 such that Rij{c*) = Rji{c*) (where Rij{c) 
is of the form Uf)) is called as a detailed balance equilibrium solution. 


2.2 The stoichiometric subspace 


Similarly to the subspace A, we define the stoichiometric subspace H as: 

^ + • • • + + • • • + 


where: 

Ha = spanj^i - \ i £ Cx} for A = 1, • • • , £ . (11) 


In what follows it will be more convenient to collect the elements from each of the sets {Vi - 
yj^ I i £ Cx] and their union, column-wise in matrices Sx £ and S £ 

respectively, so that: 


5=[5i ••• 5 a ••• 5,]. 


( 12 ) 


Let s = dim(H), which eventually coincides with the rank of 5, then it follows from the rank- 
nullity theorem that the dimension of the kernel (null space) of 5 will be: 


6 = n — £ — s. 


(13) 


This number is known in CRNT as the deficiency of the network. In a similar way, we can define 
the deficiency of each linkage class as the dimension of the kernel of Sx so that 6x = Nx — 1 — sx, 
where sx = dim(HA). Since s < it is not difficult to conclude that linkage class and 

network deficiencies relate as: 

(I>^5a. (14) 

A 
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Let {g^ I r = 1,..., J} be a basis for the kernel of S, and express each vector g'’ G ^ in terms 
of i sub-vectors g^ G (one per linkage class), so that: 

(g^)^ = [(gO^ ••• ••• (g^'^ ], for r = I,--- ,(5. (15) 

Using the above description, equation Sg^ = 0 can be re-written as: 

^5Ag); = 0 for r = l,...,6. (16) 

A 

We will be particularly interested in solutions of system ([5|) on the convex region resulting from 
the intersection of the non-negative (respectively positive) orthant in the concentration space 
and a certain linear variety associated to the stoichiometric subspace H, the result known in 
CRNT as a stoichiometric (respectively, positive stoichiometric) compatibility class. Given a 
reference concentration vector cq, the stoichiometric compatibility class can be formally defined 
as: 

0(co) = {cgM”* I c>0,P'^(c-co) = 0}, (17) 

where P G ig a full rank matrix whose columns span the orthogonal complement 

H"*-. The corresponding positive stoichiometric compatibility class can expressed as r2'*“(co) = 
n(co) n In passing, let us define the function cr : —)• as <t(c) = P^c. Such 

function is constant along trajectories dH), since by combining d?]) and ([ 2 ) we have that: 

c = XUZ ( 18 ) 

A ieCx 

and the columns of P are orthogonal to S, hence & = P^c = 0. In other words, cr is an invariant 
of motion for system dl])- From this observation it is not difficult to conclude that any trajectory 
that starts in a compatibility class ll(co) will remain there. 

2.3 Some examples of chemical reaction networks 

A reversible chemical reaction network 

Let us consider a reaction network involving m = 6 molecular species we label as {Mi,..., Mg}, 
each of them constituted by a combination of three types of functional groups (or atoms) we 
denote as A, B and C. The (reversible) chemical reaction steps that take place are: 

A 2 B + C^ AC+ AB 

AB + 2C^ AC 2 B (19) 

AC 2 B ^AC + CB 

Molecular species and functional groups are related as follows: Mi = A 2 B, M 2 = AC, M 3 = AB, 
M 4 = C, Mg = AC 2 B, Mg = CB. The network consists of n = 5 complexes: 

{Ci,C2,C3,C4,C5} = {Ml -|- M4, M2 + M3, M3 -|- 2M4, Mg, M2 -I- Mg}. 


10 


Making use of the formal description previously discussed, the sets X* that indicate which com¬ 
plexes are reached from complex i become, for this example: 


Xi = {2} X2 = {1} X3 = {4} X4 = {3,5} X5 = {4}. (20) 


The corresponding stoichiometric vectors yi associated to each complex are written as columns 
in the molecularity matrix Y: 


■ 1 0 0 0 0 

0 10 0 1 
0 110 0 
1 0 2 0 0 

0 0 0 1 0 

0 0 0 0 1 


( 21 ) 


The graph representation is depicted in Figure [TJ and comprises two linkage classes £i = {1, 2} 
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Figure 1: Graph representation for the reaction network described by reversible steps (m. 

and C 2 = {3,4,5}. The net reaction fluxes ([6]) around each complex are: 

= k2lll^2 - kl2’ll^l 
4>2{'lp) = kullJl - k2lll^2 

(()3 (■?/>) = A:43V^4 - A:34V^3 (22) 

= k^ilps -I- /C54V^5 - (A;43 -k /c45)V’4 

= kibi’A - 

Note that if ji = 1 and j 2 = 3 are chosen as reference complexes, by relation ([8|), the fluxes 
associated to the reference become: 


= -4’2W and hW = -(<(' 4 ( 1 /’) + <(>5(i/’))- 


The image of coincides with subspace A (Eqn Q) with Ai and A 2 of the form: 

Ai = span {(£2 - ei)} ^ 23 ) 

A2 = span{(£4 - £3), (£5 - £3)} ■ ^ 

Matrices S \, employed in Section 12.21 to define (column-wise) the corresponding stoichiometric 
subspaces m, are of the form: 


■ -1 ■ 


0 

0 ■ 

1 


0 

1 

1 

II 

Co 

-1 

-1 

-1 

-2 

-2 

0 


1 

0 

1- 

0 

1_ 


1 

0 

1 
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The dimension of the stoichiometric subspace, which coincides with the rank of matrix S = 
[ 5i 52 ], is s = 3. Hence, network deficiency is <5 = 5 — 2 — 3 = 0 what means that no vector 
other than the zero vector exists such that = 0. 


The numbers of atoms (or functional groups) A-C remain constant, provided that reactions 
take place on a closed domain (i.e. no mass exchanges with the environment occur). Let 
([H]o, [H]o, [CJo) be total concentrations for A, B and C on the closed and homogeneous domain. 
Mole-number balances result in the following set of linear relations: 


[HJo = 2[A2B] + [HC] + [AB] + [AC 2 B] 
[H]o = [A 2 B] + [AB] + [AC 2 B] + [CB] 
[C]o = [AC] + [C] + 2[AC2B] + [CB] 


where brackets indicate chemical species concentrations. Previous relations can be written in 
matrix form as follows: 


—co) = 0 with P = 


2 110 10 
10 10 11 
0 10 12 1 


where c is the vector of chemical species concentrations and is the (constant) concentration 
of functional groups/atoms. It must be noted that the above expression is employed in (fT7)l to 
characterize the set of compatibility classes. As discussed in Section 12.21 because rank(P) = 
m — s = 3 (full rank), the columns of P define a basis for the orthogonal complement of the 
stoichiometric subspace. 


An irreversible network 


Let us consider the following set of irreversible reactions: 


5 + 1^2/ 
S 


(24) 


This network can be interpreted as an extension of the SIR epidemic model [38j which describes 
the effect of a disease on a large population. Individuals on the population are classihed either 
as those susceptible to the disease (5), infected (/) or those recovered from the disease (R). In 
this extension, individuals under recovering may evolve either to those susceptible to the disease 
or directly infected again. In the CRNT formalism, the network comprises three species, with 
concentrations [5], [/] and [R], and 5 complexes, numbered as: 


{Ci,C2,C3,C4,C5} = {2/,5 + /,5,R,/}. 


Graph structure is depicted in Figure [2j Molecularity matrix Y for this network reads: 


Y = 


0 110 0 
2 10 0 1 
0 0 0 1 0 


(25) 
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Choosing as reference complexes ji = 1 and = 3, the S\ matrices become: 



1 ■ 


■ -1 

-1 ■ 

51 = 

-1 

, ^2 = 

0 

1 


0 


1 

0 


Because both matrices are full rank, the dimension of Hi and H 2 is si = 1 and S 2 = 2, 
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Figure 2: Graph representation for the two linkage class, irreversible reaction network (|24p . 

respectively. Thus (5i = — 1 — si = 0 and 62 = N 2 — 1 — S 2 = 0- However, matrix 

5 = ( 5i ^2 ) is rank deficient since vector 5i is parallel to the vector in the second column 
of S 2 - Consequently s = 2, and 5 = 5 — 2 — 2 = 1, verifying inequality (fT^ . A basis for the 
kernel of S is given by the vector = ( 1 0 1 )^. 

For this network, the basis that spans the orthogonal complement of the stoichiometric subspace 
corresponds to P = ( 1 1 1 )^. Each compatibility class is given by H(co) (see (fT7|) b the 

region of non-negative concentrations (c > 0) satisfying: 

[5] + [/] + [R] = P^co, 

with Co > 0 being a constant vector. A representation of a compatibility class for the reaction 
network considered is presented in Figure [3l with c = {[S], [/], [P])^. 

3 Linkage classes, Graphs and Compartmental Matrices 


The nature of equilibrium solutions in chemical reaction networks (e.g. positivity, instability of 
equilibrium points, or the possibility of multiple equilibria) is at a large extent determined by 
the graph structure of each linkage class and the properties of some matrices associated to it, 
that belong to the class of compartmental matrices |26| . 

In this section, we describe such matrices and discuss their properties, with emphasis on invert- 
ibility and non-negativity of their inverses. The main results, summarized in Lemma (3.11 will be 
extensively employed in the sequel. For the sake of completeness, we introduce a derivation from 
scratch, while establishing connections with known facts in the field of positive linear systems 
and non-negative matrices. 
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[ 5 ] 



Figure 3: Representation of a compatibility class for the network example (I24|) . 


To that purpose, we introduce a graph related to the graph description of a linkage class, and 


a matrix associated to it. The properties of this matrix will be studied by constructing an 
auxiliary linear dynamic system and examining the corresponding equilibrium. 

A directed graph Q = {{V[jvE),S} is constructed by a set of vertices containing a distinguished 
vertex ve and a set £ of edges. The first set of vertices V = {ui, U 2 ,..., Un}, with indexes 
L = {1 ,... ,n} will be referred to as nodes, while the remaining vertex ve will represent the 
‘environment’. To any directed edge Vi —>■ Vj for € £ and i ^ j \n Q (i.e. {vi,Vj) € £) there 
corresponds a scalar weight Vij > 0. In addition, for every i £ C, such that {vi^ve) € £ or 
{vE,Vi) £ £, we associate scalar weights 6* > 0 and a, > 0, respectively. Such weights will be 
collected as entries in (non-negative) vectors a, b € M”', with a* = 0 (respectively, bi = 0) if there 
is no directed edge ve ^ Vi (respectively, Vi ^ ve )■ As in the description of linkage classes 
(Section [2|) we say that two nodes are strongly linked if they can be reached from each other by 
directed paths. The maximal set of strongly linked nodes (not passing through the environment) 
in £ from which there are no outgoing edges to other nodes, defines a strong terminal set, we 
will refer to as £p. 

Since in this work we are interested in linkage classes with just one strong terminal linkage class, 
the graphs we consider will only contain one strong terminal set. The set of non-terminal nodes 
is defined as £q = £ \ £p. Some examples of directed graphs are illustrated in Figure [H We 
associate to the graph Q a state vector z £ M”, where each component Zi corresponds to a node. 
For each node i = 1,... ,n, we define an internal net flux (pi : M"" —>■ M and a net exchange with 
the environment (/jj : M ^ R. Each flux is of the form: 



(26) 
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(a) 


(b) 



Figure 4: Some typical examples of directed graphs. For both cases, the terminal set is £p = 
{1, 2} whereas the set of non-terminal nodes is /Iq = {3,4}. (a) Exchange with the environment 
takes place in the strong terminal set through input 02 (second coordinate of vector a) and 
output 61 (first coordinate of vector b). Node Env represents the environment, (b) The input 
from the environment enters the system through a non-terminal node (coordinate 04 of vector 
a) while the output to the environment leaves from a terminal node (coordinate 61 of vector b). 


where as in Section [2l the indexes of the nodes that are directly reached from node i are grouped 
in a set Ij, and {j\i € Xj} refers to all nodes j with edges directed to i. In addition, the net 
exchange with the environment is expressed as (pi{zi) = ai — biZi. From (|26l) . and similarly to 
expression ([5]), it is straightforward to see that: 

= (27) 

iec 

We consider that for each node i, the state Zi will evolve in time as a function of the corresponding 
internal and external net fluxes, so that Zi = (pi{z) -I- ipi(zi). Combining this expression with 
(I2Z1), the dynamics of the system can be described as: 

n 

z = ^4>i{z){ei-ei) - Bz + a., (28) 

i=l 

where B = I^(b), a diagonal matrix with the components of b in the diagonal. Equation (I28p 
can be re-written in the alternative form: 


z = VFz -I- a (29) 

where matrix W = V^—B, and V G is the matrix that contains as off-diagonal components 
the coefficients in (f26]l , with Vij = 0 if there is no directed edge Vi —>■ Vj. The diagonal elements 
for V are of the form: 

n 

Vu = - 

j = 1 
j + i 
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Both matrices and W are compartmental [26] , □ and belong to the class of Metzler matrices 
|6]. It is known from e.g. [7], that the eigenvalues of a compartmental matrix are either zero 
or they have negative real parts. Such conclusion can be also reached from the structure of 
compartmental matrices by applying the Gershgorin disc theorem m- 


Definition 3.1 We say that the matrix W G associated to the system is C-Metzler 

if its entries are of the form: 


Wij >0 for i j 

Wii = -{bi + with bi>0 and Wu < 0, 

with at least one positive bi associated to the strong terminal set (i.e. i € Cp). 


(30) 


Proposition 3.1 Consider system i29\) with a > 0 and nonnegative initial conditions z(0) > 0. 
Then z{t) > 0 for every t > 0. § 


Proof: In order to prove the statement all we need is to show that the flow associated to the 
differential system on the boundary of the positive orthant is either aligned to the boundary 
or oriented to the interior of the orthant. Before we compute the flow, let us define the set 
Hk = {z > 0 I = 0} which characterizes the A:-th facet of the positive orthant. The inner 
product between the flow induced by (l28l) (equivalently (12^ 1 on any element z G ff/e, and the 
unit vector orthogonal to H]^ takes the form: 

ejz = 4>k{z) - bkZk + ak = ^ VjkZj + ak>0 (31) 

lilfcGX,} 

where the equivalence at the right hand side holds since Zfc = 0 in Hk- Thus, at the boundary 
of Hk, the flow associated to the differential system will be either aligned to the boundary 
{ejz = 0) or oriented to the interior of the positive orthant (i.e. ejz > 0). Repeating the 
argument for all values of k completes the proof. □ 

Next, we will study the equilibrium of system (|29l) that results from some constant non-negative 
vectors a and b. To that purpose, let Ci C Tq be the set containing those non-terminal nodes 
that are directly linked to any node in the strong terminal set (the different partitions are 
illustrated in Figure [5|). Let us introduce functions <Tp(z) = uf^z, <Tq(z) = cjqZ and ct(z) = 
(Tp(z) -|- (Tq(z), where: 

CJp = 'y ^ £{, UJq = y ^ £{. (32) 

ie£p ie£q 

^ A matrix A G is compartmental if: (i) Aij > 0, for i,j = i j, (ii) < 0, for 

j = 

®This result is actually a special case of Theorem 2 in [20] (page 14). 
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Figure 5: Example of a graph with a, b = 0, partitioned into a terminal set £p and nonterminal 
sets £q 2 . In this case set comprises nodes 2 and 4. 


The time derivative of cjq along system (1281) is of the form: 



(33) 

iG/Zq j&Ii 



(34) 

ieCp ieCi jeXi 



Proposition 3.2 Consider system \29j) with a = 0 and W C-Metzler fDefinition \8.1\} . Then, 
for any z(0) > 0, z* = 0 is the only equilibrium solution of (2^. 

Proof: Since z(0) > 0, by Proposition 13.II we have that z{t) will remain in the positive orthant, 
cjp(z) and a^fz) will be non-negative, and cr(z) = crq(z) -|- (Tp(z) > 0. 

For every node i G ££ it is clear that z* = 0 is the only possible equilibrium. Otherwise, from 
(|33l) . it follows that dq < 0 for all times, what would make crq(z) to become negative. In addition, 
note that for z* = 0 to be an equilibrium point for node i & Ci, requires the equilibrium states 
for all nodes j linked to i (thus j G Tq) to be zero as well, otherwise Zi = <^i(z) -I- ipi{zi) > 0. 
Repeating the argument upstream to the nodes linked to j, we conclude that zero is the only 
possible equilibrium for every node in £q. 

Since W is C-Metzler (i.e. there exists at least one index i ^ Cp for which bi > 0), zero is 
also the only possible equilibrium for every node in the terminal set. Suppose, on the contrary, 
that there exists a positive equilibrium z* for some i G Cp. In the limit, expression (IMl) would 
become: 

lim dp = — > biZ* < 0, 

ie/:p 

but crp(z) cannot become negative, thus z\ = 0. Finally, since z\ = 0 is an equilibrium point 
for i G £p such that hi > 0, the equilibrium states for all nodes j linked to i must be zero as 
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well, otherwise ij = </>i(z) + ipi{zi) > 0. repeating the argument for every j G Cp we have that 
limt^oo z(i) = z* = 0. □ 

Remark 3.1 Using the terminology of j2dj. V'^ defines a compartmental system with one trap, 
where the notion of trap is equivalent to the definition of strong terminal set used in this paper. 
Therefore, zero is an eigenvalue of V'^ with multiplicity 1. W is also a compartmental matrix, 
and due to the choice of bi, the compartmental system corresponding to W contains no traps. 
Therefore, W is of full rank and thus, assuming a = 0, the only equilibrium of (flPj) is z* = 0. 

Proposition 3.3 Consider system \2fKl . with matrix W being C-Metzler fDefi,nition \S.l\} . Equi¬ 
librium z* will be non-negative and globally asymptotically stable. Moreover, let some entry k of 
vector a positive (i.e. > Oj. Then, for all nodes j ^ C reached from node k by directed paths, 

we have that z* > 0. 


Proof: First note that because matrix W in (I29p is compartmental, its eigenvalues are either zero 
or they have negative real parts. Thus, to show that the equilibrium is globally asymptotically 
stable it only remains to prove that no zero eigenvalue exists. Actually, this is the case, since 
from Proposition 13.21 the only equilibrium solution for which VFz* = 0 is z* = 0. Since all 
eigenvalues have negative real part, W is invertible and the equilibrium z* = —W~^a is globally 
asymptotically stable. 


In proving the second part of the statement, we have that since > 0, the equilibrium in node 
k must be: 

, EfaltsJ,) CfiZ,' + at 
‘ + Eia. Vtj 

and so is the case for all j reached from k, so that z* > 0, what completes the proof. □ 


Remark 3.2 The full rank property ofW and the uniqueness of the equilibrium zf" = 0 for a = 0 
are equivalent, since this latter means that the dimension of the kernel of W is zero. From this, 
it follows that W cannot have zero eigenvalues, and from the compartmental property of W we 
obtain that all the real parts of its eigenvalues are negative. This means that W is a stability 
matrix. Since a can be considered as a bounded input, z* = W~^{—a) is a unique asymptotically 
stable equilibrium point of (2^. Since —W is a full-rank M-matrix, it is inverse-nonnegative 
all entries ofW^ are non-positive. This implies that zf" is a nonnegative vector. 


Lemma 3.1 Any C-Metzler matrix W is non-singular and its inverse W~^ non-positive. Let 
its associated graph he numbered so that the first p nodes are in £p (thus, the remaining q = n — p 
are in £qj. Then N = W~^ can he partitioned as: 


N = 


-1 


0 

\ 


(35) 
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with 0 G the zero matrix, N-p G RP^p and A^pq G RP^p strictly negative matrices, and 

Nq G non-positive. Moreover, for each node j G Cq, every node i not reached from j by 

directed paths will correspond with an entry Nij = 0. 


Proof: That W (being C-Metzler) is non-singular and therefore invertible has been shown in the 
proof of Proposition I, I., 'll (see also Remark 13.21) . In addition, note that the equilibrium solution 
can be written as z* = —Na, and choose a = sj. For each j G C, the corresponding equilibrium 
will be z* = —(N)j, where {N)j represents the j-ih. column of matrix N . 

According to Proposition [331 the equilibrium state for all nodes i G T reached from j by directed 
paths must be positive, and thus the corresponding entries Nij must be negative. This holds 
for all i,j G Cp. Thus, using the order given in the statement of the Lemma, we can conclude 
that Np is strictly negative. In addition, every i G Tp can be reached from j G Cq what makes 
A^pq strictly negative as well. The zero matrix 0 appears since none of the nodes z G Tq can be 
reached from nodes j G Tp. Finally, zero entries in A^q will correspond to those i not reachable 
from j, with i,j G Cq. □ 

Lemma [3.II will be invoked along the sequel on a linkage class basis to study positive equilibrium 
solutions of chemical reaction networks. The possibility of positive equilibria is at a large extent 
connected with Proposition 4.1 in Lecture 4 of Feinberg [22] (see also proofs in [3l| and IZSj) on 
the structure of the kernel of defined in Q. In this regard, it is noted that the arguments 
behind Proposition 13.31 and Lemma 13.11 can serve as a basis to prove the result in Feinberg 
Lectures. 

Example: The graphs in Figure 0 ] depict two different scenarios. In Figured^, input 02 enters a 
node in the strong terminal set thus forcing the states z* in that set to be strictly positive, while 
leaving those that belong to the non-terminal set to reach zero. Figure [Hd, describes the case in 
which the input enters a non-terminal node that communicates with the remaining nodes of the 
graph (non-terminal and terminal), forcing the system to reach a strictly positive equilibrium 
state z* > 0 . 


Since 61 > 0 leaves a node from the strong terminal set, the associated matrix W is C-Metzler 
fDefinition 13.ip . The sign patterns for the first two columns {N)i and {N )2 of its corresponding 
inverse N are of the form: ( — — 0 0 )^. For the 3'”'^ and 4*^ columns, the sign patterns are 

( — — — — )^, since every node can be reached from either node 3 or 4. The sign pattern 

of the negative inverse will then be: 


i-N) 


+ + + + 
+ + + + 
0 0 + + 
0 0 + + 


A 
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4 A Canonical Representation of the Equilibrium Set 


The time evolution of the concentration vector (|18p can be re-written in terms of the Sx matrices 
associated to the stoichiometric subspace H, already discussed in Section [221 so that: 

c = '^Sx(l)x{i^j^{c),'tpx{c)), (36) 

A 

where : ^>0 I^>0) corresponds to the monomial associated to the reference complex, and 

: IK>o vector function that includes as coordinate functions the monomials 

associated to the remaining complexes. Finally, vector function (px ■ lK>o contains 

the corresponding fluxes. Element-wise, monomials relate to fluxes by expression ([6]), that in 
matrix form can be written as: 


(37) 


where the fluxes and monomials are ordered such that for each linkage class, the first elements 
correspond to the reference complex ja (i.e. ^jy(c),?/)j^(c)). Matrix Mx G is a compart- 

mental matrix that has as off-diagonal entries the corresponding reaction constants. An explicit 
description of its structure can be given as follows: 


Let Cx{i) denote the i-th element in the set £a for i = 1 ,..., Nx- Without loss of generality, we 
can assume that the first element of Cx is the index of the reference complex, i.e. /1 a(1 ) = jx- 
In addition, let kjr^(i)^Cxij) denote the reaction rate coefficient associated to a possible reaction 
step from complex Cx{i) to Cx{j), being 0 if such reaction step does not exist. Then: 

{Mx)ij = kcx{j),Cx{i) hi = 1,---,^A, (38) 

Nx 

{Mx)ii = - ^ kc^^i^^CxH) fo'^ f = l,...,A'A. (39) 

1 = 1 , 


Note that because of dEI), we have that: 

</>ja(^ja>'0a) + l%^_i(l)x{'>Pjx,'fPx) = 0, 


(40) 


so that Mx is a Kirchoff (i.e. 
as: 


column conservation) matrix, which for convenience we re-write 


Mx 


(1?^x-1^a) 

bl 1 

^A 

Ex \ 


(41) 


with aA,bA G Ex G and = — 1^^_^Ea. By construction, the off- 

diagonal elements of the first column and row in Mx correspond to the rate coefficients for 
reaction steps leaving and entering, respectively, the reference complex jx- Such reference has 
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been chosen to be in the terminal linkage class, what ensures to be non-zero, since at least 
one reaction step is directed to the reference complex. The off-diagonal elements of matrix Ex 
collect the remaining rate coefficients. 

Proposition 4.1 Ex in expression eip is C-Metzler, therefore invertible, and its inverse non¬ 
positive. 


Proof: First, we note that Ex complies with Definition 13.11 (C-Metzler matrices). This is so 
because it is associated to a directed graph which coincides with the linkage class (the environ¬ 
ment corresponds with the reference complex). By construction, the off-diagonal elements of Ex 
are either zero or positive. In addition, since Ex, for each diagonal element we have 

that: 

{Ex)ii = ~[(bA)i + . fEx)ji], 

where at least one component of b;^ associated to the reference complex (linked to the strong 
terminal set) is positive. Thus, Ex is C-Metzler and the result then follows by applying Lemma 

EH □ 


Inspection of Eqn (l36]l suggests that apart from complex balanced equilibrium solutions (i.e. 
those satisfying Definition EH), non-zero ffux combinations can lead to equilibrium, if the cor¬ 
responding flux vector cf>{rp) = [ (pj , '0i) • • • (p'x > '*/’a) • • • (pj{'f’je > '^i) in the 

kernel of S (I12p . In other words, if the network has non-zero deficiency, vector can be 

written as a linear combination of the set of vectors {g^ | r = 1,..., <5} that dehne a basis for 
the kernel of S, so that: 

(Pi'fP) = (42) 

r 

for some given scalars Making use of (|15l) in the above summation to express each element 
r of the basis in terms of the sub-vectors g^, for A = 1, • • • ,i, the flux vector for each linkage 
class can be written as: 

0A(V'jx>^A) = (43) 

r 

Note that such fluxes lead in fact to an equilibrium solution. This can be verified by substituting 
(in|l into (PS]1 and re-grouping summations, so that: 

c = '^Sx'^i^rgl = '^i^r'^Sxgl, (44) 

X r r X 

where the right hand side is zero because of ()16p . 


Example: Let us consider the network depicted in Figure EJ Since £ 


1, there is just one 


matrix Mi (1411) that consists of the following sub-matrix components: 


-(k21 + k23) 

k32 

0 

0 ■ 


0 ■ 


k2i 

k23 

— k32 

kis 

0 


0 

, bi = 

0 

0 

0 

”(^43 + ^45) 

^54 

, 3-1 

0 

0 

0 

0 

^45 

— ^54 


ki5 


0 


(45) 
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(a) 


(b) 




Figure 6: A one linkage class weakly reversible network consisting of 4 chemical species 
{A, B,C, D}. (a) The graph of complexes with explicit indication of the chemical species, (b) 
The graph with numbered complexes, including the reference complex in the square box and 
the graph associated to the C-Metzler matrix with the non-zero coordinates of vectors ai and 
bi, represented in the diagram as (ai)4 and (bi)i, respectively. 


The Y and S matrices described in Section [2] are, respectively: 


■ 1 

0 

0 

2 

0 ■ 


■ -1 

-1 

1 

-1 ■ 

0 

0 

0 

0 

1 

and S = 

0 

0 

0 

1 

1 

2 

0 

0 

0 

1 

-1 

-1 

-1 

_ 0 

0 

1 

0 

0 _ 


0 

1 

0 

0 _ 


(46) 


Since rank{S) = 3, the network has deficiency <5 = 1. Thus, equilibrium solutions would 
correspond to fluxes satisfying relation (1431) with g^ = g( = (l 0 1 O) defining a one¬ 
dimensional subspace. A 


In computing the set of positive equilibrium solutions, we make use of (l37)l and m for every 
linkage class, to express fluxes cj>x as: 

(47) 

Combining the right hand sides of (I43p and (|47D we get: 

Ex-il^x + ^rg\- (48) 

r 

Let u = [i/i ... and define a vector % € on the unit sphere as % = Since any 

positive equilibrium solution requires 'ipj^ > 0 and tpx > 0, we can multiply both sides of Eqn 
(1481) by I/'jAj'a obtain, after some re-arrangement: 

Exh + aA = xxGxX, with fA = {l/'tpj^)'il7x- (49) 

In the above expression, Ga G is a matrix of the form Ga = [gA •••§)(••• SaIj and 

xx = a scalar variable. It must be noted that because = xxil’i. for every A = 1,... i, 

Y3\ 

variables xx are not independent but related to each other through the following equalities: 

Xi'i/’ji = • • • = XA'i/’iA = • • • = Xt,i}j^. (50) 

Since Ex is invertible (Proposition 14.1|) . we can solve (I49p to get a vector function Ta : M x —)■ 

]^(AfA-i) Qf form: 

fA(a:^A;A:) = + ®AhA(A:), (51) 
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where and hA(x) = ^gA(x)) with gA(x) = G\X- Vector fluxes in (H7|) can be 

expressed as a function of (fCT) for each linkage class so that; 

4>\ii’jx,xx-,x) = V'ix [aA + -EAfA(a:A;x)] • (52) 

In this way, for each x ia the unit sphere and x = (xi, • • • ,x\, ■ ■ ■ ,xi)^ , constrained by (l50]l 
to be either zero or to belong to the interior of the positive orthant Kio: the right hand side of 
(jlibp vanishes. This follows since for each x on the unit sphere, the corresponding fluxes (|52p 
become: 

(t}x{'ipjx,xx;x) = x\ipj^GxX- 

Substituting the above expressions in (pbl) and making use of (fSOl) we get: 

where the right hand side is zero because from (uni) we have that: 

^aGa = 0. 

A 

Note that for a vector x' = ~X we have that x € {0} U R'<o. SO it is enough to study vector 
functions (f5T]l for x € {0} U Ml.0 U M^o and X ia the set: 

U = {xGM'I llxll = 1, efx>0}. (54) 


Definition 4.1 (Family of Solntion^) Let Bq C be the domain {0} U M^q U K^g, and 
fr, : Bo X U —)> a vector function: 

f^(x;x:) = [ ff (xi;x) ••• i'[{xx;x) ■■■ ^iixe-,x)], (55) 

with X = (xi,--- ,X£)^ and fA(a;A;x) for X = 1, - ■ ■ ,i as in We will refer to f^(x;y;) as 

the family of solutions. 

In order to comply with positive equilibrium solutions in the concentration space, each vector 
function fA(xA;x) in (f5Tp must be strictly positive, and related to the reaction monomials by 
the expression: 

fA = exp ^In :^V’A^ • (56) 

In what follows, and at the risk of some abuse of notation, when referring to a given linkage 
class, we will drop subscript A, and re-write m as: 

f(x;x) = r-bxh(x:). (57) 

closely related notion has been proposed in m 
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If in addition, the discussion concerns a particular vector % G U, the following simplified ex¬ 
pression for Eqn ([Fr|) will be employed: 


f(x) = f* -|- xh. (58) 

Lemma 4.1 If the linkage class is weakly reversible, f* in is strictly positive. If the linkage 
class is irreversible, f* will be non-negative with zero entries that correspond to the non-terminal 
complexes. 

Proof: If the linkage class is weakly reversible, every complex in the linkage class can be reached 
from the reference (equivalently, from the complexes associated to positive entries a in (1411) 1. 
Thus from Lemma [TTI f* = —E“^a > 0. 

If the linkage class is irreversible, let a^ = [ap aq] where sub-indexes p and q denote the 
terminal and non-terminal complexes of the linkage class. Since the reference belongs to the 
strong terminal linkage class, positive components of vector a only enter terminal complexes so 
that aq = Oq. Using the inverse (|35l) from Lemma l3.11 P = —Na. can be written as: 


r p 1 

p 


-A^pap 

p 

L J 


. Oq . 


where ^ is strictly positive because, by Lemma IXTl Aip is strictly negative. Thus P > 0, with 
zero entries f^ that correspond to the non-terminal complexes. □ 

Next we present some conditions that ensure positivity of vector functions (|5ip and the corre¬ 
sponding family of solutions (Definition 14.ip . 

4.1 Positivity conditions for the family of solutions 

At this point, it should be clear that positivity of the family of solutions in Definition 14.11 is a 
necessary condition for positive equilibrium in the concentration space. Next, we discuss how 
such condition relates to the structure of the network and give some indications on how to 
construct the domain where (|55p remains positive. 

Proposition 4.2 If a given linkage class is weakly reversible, then for every y G U there exists 
an interval X{x) C M (which includes the zero) where the vector function remains positive. 
If the linkage class is irreversible, and there exists a positive vector function (5^ on X(x), such 
interval cannot contain the zero. 


Proof: If the linkage class is weakly reversible, by Lemma [4. II we have that P > 0, so the values 
of the scalar x for which f(x) in (1581) remains strictly positive will depend on the signs of the 
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entries in h. For every such entry i, define pi = hi/f* and introduce two index sets and Z 
so that: 

i G if hi > 0, thus Pi > 0 
i G Z~, if hi < 0, thus pi < 0. 

It is straightforward to see that f(x) will be strictly positive for every x in the open interval: 


(59) 


X = (L with L = maXjg2;+{~l/Pi} and, L'*'= minjgj-{—l/pi}. 


(60) 


Note that L~ = —oo (respectively, L+ = +oo) provided that = 0 (respectively, Z~ = 0). In 
any case, because f* > 0, the interval X includes the zero. If the linkage class is irreversible, by 
using Lemma l4.11 f(x) in (1581) can be written as: 


fp(x) = + xhp 

fq(x) = Oq + xhq. 


(61) 


where sub-indexes p and q denote the terminal and non-terminal nodes of the linkage class. Since 
^ is strictly positive, there exists some domain Xp that includes the zero, for which fp(x) > 0. 
Let Xp = Xp UXp U{0}, where X” and Xp are the intervals containing the negative and positive 
values, respectively. It is then straightforward to see from (fST]l that in order for fq(x) > 0, hq 
must have a definite sign (i.e. all components either positive or negative). If this is the case, 
i.e. ifhq > 0 (respectively, < 0), we can always hnd some x G X+ (respectively, x G Xp), so 
that f(x) > 0. Otherwise, no positive solution exists. Because X+ (respectively, X" ) does not 
contain the zero, if the linkage class is irreversible, the interval X(x) does not contain the zero. 
□ 

Proposition 4.3 Letfri : BqxU ^ (withOo = {0}UM^qUM^qJ be the family of solutions 

as given in Definition \4.1\ In addition, for each A = 1, • • • and x G U, let I^x{x) C M 6e the 
interval such that Ixixx-,x) > 0 for every xx G Xa(x), andXr,{x) = Xi(x) x • • • Xa(x) x • • • X£(x) 
an open i-dimensional domain. Then: 

If the network is weakly reversible, for every X G U there exists a domain ID)^(x) = X^(x) H Dq 
that contains the zero, such that fr,(x;x) > 0 for every x G Bj,(x). 

If the network is irreversible and the domain B^(x) = X^(x)nBo is non-empty, then frj(x; x) > 0 
for every x G Bj,(x). Such domain does not contain the zero. 

Proof: If the network is weakly reversible. Proposition 14.21 ensures that for every A = 1, • • • 
and X € U, fA(iCA;x) > 0 for every xx G Xa(x)j where the interval contains the zero. In 
consequence, a non-empty domain Bj,(x) = X^(x) H Bq exists, what proves the first assertion. 


If the network is irreversible and B^(x) is non-empty, the second assertion holds. Because at 
least one linkage class is irreversible, it follows from Proposition 14.21 that X^(x), and therefore 
B^(x), does not contain the zero. Finally, note that for irreversible networks it may well happen 
that Bj,(x), as defined above, is empty so no positive family of solutions exists. □ 
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4.2 The set of feasible (equilibrium) solutions 


Not all positive elements (vectors) that are part of the family of solutions (|55p will necessarily 
comply with condition (l56|) but only a particular subset we will refer to as the set of feasible 
solutions, that we formally define next. Using ([2]) for every complex i € C\, we have that 
ln('0j(c)/'!/)j^(c)) = {ui — yj^)'^\nc. Hence, the logarithm at the right hand side of (f56]l can be 
expressed as: 

In , \ x ^a(c) = Sj Inc, (62) 

Vjx\^) 

so that Inf^P^) = for every A = 1, and ^ € M”* = Inc). 

Definition 4.2 (The set of feasible solutions) For a given x € U, let us assume that there 
exists a non-empty domain lD>jj(x) C Bq such that for every x G Bj,, f, 7 (x; x) > 0. We say that 
fr,(x;x) > 0 is a feasible solution if: 

lnf^(x;x) G Im(S''^). (63) 

The set of vectors fr,(x;x)j with x G U <ind x G Bj,(x) that satisfy ^6^, constitutes the set of 
feasible solutions. 

Note that (1651) implies that there exists some ^ G such that lnfj,(x; x) = S'^^. In this way, 
each element of the set of feasible solutions relates to a set of equilibrium concentrations of the 
form c = exp(^) for system (|36p . The following result gives some conditions to identify the set 
of feasible equilibrium solutions; 

Lemma 4.2 (Feasibility conditions) Every element f^(x;x) of the set of feasible solutions 
satisfies that: 

(g^)^lnf^(x;x) = 0, for r = l,...,6. (64) 

Proof: By Definition 14.21 every element of the set of feasible solutions satisfies Eqn (|63p . This 
implies that lnfjj(x;x) is in the range of S'^, which in turn is orthogonal to the kernel of S. 
Since {g^ | r = 1,..., d} is a basis for the kernel, expressions in (1641) follow. □ 

Definition 4.3 (Feasibility function) Let F : X x U ^ M (with X C M) be defined as: 

F(a:;x) = g'^(A:)lnf(x;x), (65) 

where X G U and X(x) is the interval in which f(x;x); of the form ( [57| j, remains positive 

We make use of the above definition to present an immediate consequence of Lemma 14.21 
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Proposition 4.4 For every element f^(x;x) of the set of feasible solutions, the following rela¬ 
tion holds: 

= 0, where Fa(xa;x) = pI(x)lnfA(xA;x). (66) 

A 

Proof: Pre-multiplying each equality in (|64l) fLemma l4.2D by the corresponding coordinate Xr, 
taking the summation to 6 and expanding over linkage classes, we get: 

and ^gJ(x)lnfA(xA;x) = 0. 

r X 

The result then follows by using Definition 14.81 to re-write the above expression as in (1661) . □ 

As we will see in the next sections, feasibility functions FA(ic;x) will prove to be fundamental 
to characterize the structure of equilibrium, allowing in some instances to conclude uniqueness 
of equilibrium in each positive stoichiometric compatibility class. 


4.3 Example: feasibility for a one linkage class weakly reversible network 


(a) (b) 



Figure 7: A weakly reversible one-linkage class network with reaction constants ki 2 = ^34 = e 
and k 23 = k^i = ft. (a) The network with explicit indication of the chemical species, (b) The 
corresponding graph of complexes with the reference complex in the squared box. 


Let us consider the weakly reversible reaction network taken from |23j and presented in Figure 
[71 Matrix m for this network takes the form: 


M = 


—ki2 0 0 kii 

ki2 -k23 0 0 

0 k23 -k3i 0 

0 0 k3i —k^i 


(67) 


Substituting the reaction constants given in Figure [71 leads to the following sub-matrices in M: 


( 68 ) 



■ -/3 

0 

0 ■ 


e 


■ 0 ■ 

E = 

/3 

—e 

0 

, a = 

0 

, b = 

0 


0 

e 

-/3 . 


_ 0 _ 


. /? . 
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For this example, matrix S and a basis for its kernel, expressed as columns of matrix G read: 


- 1 

1 

to 

1 

1—1 

1 

CO 

1 

, G = 

■ 1 0 -2 ■ 

2 

3 1 

0 1 -3 


In order to compute the feasibility function (f6^ . we have that: 


r =-E-^ai= { e/^ 1 e//3)^, E-^ = 

r 1 

_! 

0 

_ 1 

0 ■ 

0 

1 

1 

1 


L /3 

/3 

4 J 


and 

h{x) = E-^ 

For parameters (5 = 2 and e = 



/ 1/2 
r= 1 
V 1/2 



1 we have that: 


j and h(x) = xi 


( \ 

Xi -1/e +X2 

\ 1//3/ 


- 1/2 \ / 0 

-1 +X2 -1 

1/2 / V 1 



(69) 


For xi = (1 0)^5 the entries of the vector function f(x; Xi) will remain all positive as long as the 
values taken by x will lie in the open interval (—1,+1). Hence, its domain X(xi) = (—1,+1). 
In the same way, for X 2 = (0 1)^, ^{X 2 ) = (“l/2,+l)- Because the logarithm is defined for 
positive values, both domains X(xi) and X(x; 2 ) coincide also with the domains for F(x; Xi) and 
F(x;X 2 )- The explicit expressions become: 

, 2(1 — x) 1 -n/ \ 1 8(1 — x) 

F(x;xi) =ln-^^^^, and F(x;X 2 ) = ln 

Feasibility functions for some vectors x iii the unit sphere are presented in Figure [HI It must 
be observed that each vector x leads to a different domain X(x). Finally note that since F is 
strictly positive, all possible domains will include the zero. A 


It must be noted that the functions depicted in the above example are monotonous decreasing in 
their respective domains. Remarkably, this will be the case for any feasibility function F(x;x), 
despite network structure or stoichiometry. Next section provides a formal proof of this fact 
which will turn out to be central in exploring the nature of equilibrium solutions. 


5 Monotonicity of feasibility functions 


Here, we study the properties of function (f65|) . presented in Definition 14.31 In particular, it will 
be shown that it is monotonous decreasing in its domain X(x) and crosses the x-axis. 

Let us consider a weakly reversible linkage class and a given vector % € U. Without loss of 
generality, assume that the N — 1 components of the vector h in (I58|) . being m of them positive. 
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Figure 8: Feasibility functions F{x] x) associated to the network presented in Figure [7] for 
different vectors x the unit sphere, (a) X = (1 0)^, (b) X = (0 1)^, (c) x = (“1 0)^; (d) 
X = (0 — 1)^. Note that since vectors in plots (a)-(c), and (b)-(d) relate as x! = ~X^ the 
corresponding functions relate as F(x;x) = —F(—x;xO- 

r zero and the remaining negative, are ordered so that: 

h-i ^ ^ h]^ ^ ^ hfYi > 0 > h-in+r+i ^ ^ > • • • ^ /i/v"—1, (70) 

hm+l = • • • = hm+r — 0. 

Note that such order can always be induced by a suitable row and column permutation in 
equations g = Eh and EE = —a. In order to simplify notation, let us re-write function (I65p for 
a fixed x a-s: 

F(x) = g'^lnf(x). (71) 

The main result on monotonicity is presented in the theorem below. 

Theorem 5.1 Let X C M and consider the function F(x) : X i—>■ M defined in F(x) is 

monotonous decreasing in the interval X, defined as in (Eg). Moreover, 

lim F(x) = +00 and lim F(a:) = —oo. (72) 

x+^L~ x~^L+ 


Proof: Function (j7ip is continuous differentiable in the the interval X since f(x) is strictly 
positive (see Proposition 14.2D . Thus, the first part of the proof reduces to computing the first 
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derivative with respect to x and studying its sign in the interval X. For every entry i of vectors 
h and f*, let us define pi = hi /as in the proof of Proposition 14.21 and re-write h as: 


h = v{r)p, 


(73) 


where vector p £ ^ includes the elements pi and D(P) represents a diagonal matrix with 

the components of P in the diagonal. Let us also re-order the pi elements so that: 


Pi > ■ ■ ■ > Pk> ■ ■ ■ >Pni> 0 > Pm+r+l > ■ > Pi > ■ > PN-1, 

Pm+1 = • • • = p^n+r — 0 - 


(74) 


Note that the number of positive, negative and zero elements must coincide with those in (|7nii . 
although not necessarily in the same order. Define functions Qi{x) : X i—)■ M as: 


Qi{x) 


Pi 

1 + xpi' 


For every k such that pk > 0 and x € X, we have that x > —{1/pk), since from (inU|l : 


(75) 


x > maXjgi+{-l/pi}(= -1/pi). 


In turns, this implies that x + (l/pk) > 0, and: 


Qk{x} 


1 

X + (1/Pfc) 


> 0 . 


Using the same argument for the negative elements, we have that x < —1/pi so that Qe{x) < 0 
for any ^ = m-|-r-|-l,...,iV — 1. In addition, for any pi > pj (both, either positive or negative), 
we have that Qi(x) > Qj{x). Consequently, from for every x G X we have that: 


Ql{x) ^ ^ Qk(.X^ ^ ^ Qmix) > 0 > Qrn+r+lix^ ^ ^ Qi{x^ ^ ^ Qn—i{x^i 

Qm+l{x') — • • • — Qrn+r^X^ — 0 . 

(76) 

Keeping the order established in (17411 . the first derivative can be written as: 


N-l 


F'(x) = 


9ihi 


(77) 


where gi is the i coordinate of vector g = Eh., and fj(x) represents the i component of vector 
function ()58p . The derivative is well defined and continuous on X, since fi(x) > 0 for every i and 
X € X. By dividing every element of the summation by f*, and using pi = hi/f*, we can re-write 
dZZl) in term of functions Qi{x) (USD as: 


N-l 

F'(x) = ^ PiQiix). 

i=l 


Let us define a matrix H € i)x(Ar i) 


(78) 


H = EV{r), 


(79) 
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which by construction is C-Metzler (Definition [ST]) , since E is C-Metzler, and the columns of E 
are scaled by a positive diagonal matrix. By means of H and Eqn (1731) . we re-write Ef^ = —a 
and g = Eh, respectively, as: 

= —a, and g = Hp. (80) 

Because H is C-Metzler, and the relations dZSl) and ([8(1]) hold, we are under the conditions of 
Lemmas lA.ll and IA.2I in Appendix In particular, the right hand side of (|78p has the same 
structure as G{x) in Lemma lA.21 Consequently, the hrst derivative is strictly negative on the 
interval X and monotonicity follows. 


In order to prove (f72]l . we note that each entry of f(x) can be expressed as fj(x) = f*(l -|- xpi), 
and re-write (1711) as: 

N-l N-1 N-1 

^ ffilnfj(x) = ^ c/jlnfi -h ^ gi^i{x), (81) 

2 = 1 2 = 1 2 = 1 

where nj(x) = ln(l -|- xpi). In addition, let us re-write the sequence of positive parameters in 
dZll), in the equivalent form: 

Pi = ■ ■ ■ = Ps > Ps+i > ■ ■ - Pk > ■ ■ ■ > Pm > 0, 

with s being an integer that denotes the hrst strict inequality in the sequence (counted starting 
from the largest element), and can take any value between 1 and m. Similarly, let us re-write 
the sequence of negative parameters in (f7^ as: 

0 > Pm+r+i > ■■■ >Pi> ■■■ > Pt-i > pt = ■■■ = PN-l, 

with t being an integer that denotes the hrst strict inequality in the sequence (counted from the 
smallest element), and can take any value between m -|- r -|- 1 and N — 1. 


The hrst term at the right hand side of ()8ip is constant, while the second term can be expanded 
as in (1A.14P (proof of Lemma IA.2D with Tii{x) instead of Qi{x). Taking into account the above 
sequences, the expansion can be written as: 

N-l 

F(x) = ^5ihi^ + F+(x) + F-(x), (82) 

2 = 1 


with: 


2=1 


2=1 


F+(x) = (ni(x) - ns+i(x)) H-h (nfc(x) - nfc+i(x)) ^5* 

m 

T n)^(x) 'y 

i=l 

F“(x) = Iirn+r+l{x) ^ Qi ^ -h (n£(x) - II^.i(x)) ^ 

2=rr2+r+l 

7V-1 

-F (nAr-i(x) - nt_i(x)) ^ 


(83) 


Af-1 


N-l 


i=i 


(84) 


i=t 
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where terms of the form nj(x) — Hj (x) such that pi = pj have been dropped from the expansion, 
as they are zero. In computing the left and right limits of F(x), these are the possible scenarios: 

1. If s = m or t = m + r + 1, expansions (f55|l or ([MD reduce to: 

m N—1 

F+(x) = nm(x) or F“ (x) = Ilm+r+i ^ Qi- 

2=1 2=m+r+l 

Bv Lemma |A.1| lExpression in (|A.6pL we have that Oi < 0. and r^+r+l 9i > 0. We 

also have that pm = Pi and pm+r+i = Pn-i- Thus, using the limits (IA.18h in Proposition IA.2I 
we obtain (172]). 


2. If s < m and t > m + r + 1, there are both, positive, zero and negative entries so the interval 
becomes X = (—I/pi, — I/pAr_i). In the limit as x+ —)• — I/pi all terms in (1821) are constant (see 
limits (|A.2ip in Proposition IA.2P except the first term associated to F'’'(x). Concerning this 
term, we have that Yl,i=i9i < 0 (i-®- strictly negative according to Proposition lA.ip so by using 
()A.19p of Proposition lA.21 we get: 

S 

lim F(x) = lim (IIi(x) — If^+iPx)) > gi = +oo. 

x+^L- x+^-(^) ^ 

' Cl 7 — 1 


In the limit as x — —1/pN-i, all terms in (]52p are constant, but the last one associated to 
F“(x). Again, using Propositions I A. ll and [A?2l we have that: 


lim F(x) 

x--^L+ 


N-l 


lim 




PJV-I 


(nAr_i(x) - nt_i(x)) 

) 


= -OO- 

i=t 


3. All entries are positive or zero, so that m + r = A^ — I. In this case, the interval X = 
(—I/pi,+oo), and (l82P reduces to: 


N-l 

F(x)= <7*111 fT+F+(x). 

2 = 1 


As in the previous case, \\mx;+^L- F(x) = +oo. On the other hand, in the limit as x —>■ +oo, by 
Proposition lA.21 all terms but the last one in F^(x) are constant. Because YllLi9i < 0 (Lemma 
lA.ip . and lima;_).+oo Ilm(x) = +oo, we have that: 


lim F(x) = lim nm(x) > Qi = —oo. 

X—>-+00 X^ + CXD ^ 

2=1 


4. All entries are negative or zero so that m = 0, the interval X = (—oo, —l/pAr_i) and (18211 
reduces to: 

N-l 

F(x)= ^< 7 *ln^+F-(x). 

i=l 
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As in case 2, \\ra.^-^i^+ F(a:) = —oo. On the other hand, in the limit as x —>■ —oo, by Proposition 
IA.21 all terms but the first one in F“(x) are constant, and because Tld=m+r+i9i > (Lemma 
lA.ll) and lima,^_oo n^+r+i = +oo, we have that: 


lim F(x) 

X—>• —OO 


N-l 


lim nm+r+i(a:) 

x—>• —OO 


^ gi = + 00 . 


2=7n+rH-l 


□ 


6 Some network classes with unique equilibrium solutions 


Chemical reaction network structure, with its associated C-Metzler matrices, influences the set 
of feasible equilibrium solutions that we compute by applying conditions in Lemma 14.21 and are 
related to the feasibility functions given in Definition 031 In this section, we exploit monotonicity 
of these functions to explore existence and uniqueness of equilibria within positive stoichiometric 
compatibility classes for weakly reversible reaction networks. 

Monotonicity will be used to identify sub-sets within the space of possible reaction rate coeffi¬ 
cients leading to complex balancing and in line with the classical works in [351 IMl I22| , to conclude 
existence and uniqueness of equilibria within positive stoichiometric compatibility classes. 

The monotonicity argument will also be employed to show existence and uniqueness of equilib¬ 
rium solutions for a class of positive deficiency networks. This might support the construction 
of an alternative proof of the well-known deficiency one theorem [23 [8] for weakly reversible 
reaction networks. 


6.1 Zero flux conditions and complex balanced equilibrinm 

Here we examine the equilibrium solutions c* for the dynamic system (|36p . that result from all 
fluxes in the network to be zero, namely = 0 for every \ = 1,. .. ,i. In exploring 

such a case (we will refer to as the zero flux condition), we first note that frj(0;x) (i.e. the 
family of solutions (Definition 14.ip evaluated at x = 0) corresponds to a zero flux condition. 
This can be shown by substituting = 0 for every A = 1,..., I" in ([52p . and and using the fact 
that f^(0;x) implies that fA(0;x) = f)(, so that Eqn (1521) becomes: 

4>\{ipjx,0-,x) = V'iA(aA T-FaIa) = 0, 

where the zero flux condition follows since = —E^^ax for every A = Let us denote 

f^(0; x) by f*, which by construction is of the form: 

(f*)^ = [ ••• (f;)^ ••• (f;)^]. (85) 
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If the network is irreversible, it follows from Proposition 14.31 that the domain Bfj(x), is either 
empty for some % € U, or if not, it does not contain the zero. Hence, f* cannot be a strictly 
positive vector, what in turns results in some species concentrations (associated to the zero 
entries of f*) to be zero. Since there is no strictly positive equilibrium vector c* complying 
with a zero flux condition, irreversible networks do not accept complex balanced equilibrium, 
according to Definition 12.11 

On the other hand, if the network is weakly reversible. Proposition 14.31 asserts that for every 
X € U, there exists a domain B, 7 (x), which contains the zero, such that fr;(x;x) is strictly 
positive, and consequently f* > 0. If in addition, f* belongs to the set of feasible solutions 
(Definition 321)) there exist strictly positive vectors c = exp(^), such that Inf* = 5^^, which are 
equilibrium solutions of system (13611 . According to Definition 12.11 those equilibria are complex 
balanced. 

We recall from Eqns (HT|) , (14^ and (l5T|) that f * depends on the reaction rate coefficients through 
E\ and a\ for A = !,...,£, which in the last instance determine feasibility, in the sense of 
Definition 14.21 For convenience, we collect the set of reaction rate coefficients of the network 
into a vector k € Rto. where p denotes the number of irreversible reaction steps in the network, 
and introduce the so-called Horn set pQ, that is formally defined as: 

'W = {fce]R^o I Inf* G Im(5^)}. (86) 

Note that the set is only meaningful for f* > 0, which as discussed above requires the network 
to be weakly reversible. The result we present next shows that the set H contains all possible 
reaction rate coefficients leading to complex balanced equilibrium solutions. 

Proposition 6.1 Any chemical reaction network with k £ H will only accept complex balanced 
equilibrium solutions. 


Proof: For any G "H, fj) is an element of the set of feasible solutions (Definition 14.21) . since 
there exist vectors ^ G M™' such that: 

lnf* = 5^^, (87) 

and as discussed above, the corresponding strictly positive vectors c = exp(^) must be complex 
balanced equilibrium solutions satisfying: 

Inf* = 5'^ Inc. (88) 

In fact, as we will prove next, these are the only possible equilibrium solutions. First, we note 
that by Proposition 14.41 for any x G U we have that: 

^Fa(0;x)=0. (89) 

A 
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Suppose that for a k ^ T-L there exists a non-zero flux condition that lead to an equilibrium 
solution. This implies that there exists at least one G U and x* G ror 7 (x*) with x* 7 ^ 0 such 
that fr,(x*;x*) is feasible, and therefore relation (l 66 l) applies, so that: 

Efa(4;x*) = o, (90) 

A 

Since the network is weakly reversible, by Proposition 14.31 we have that the domain 0^(x*) = 
Id Bo is non-empty, contains the zero, and can be partitioned as: 

b„(x*) = b-(x*)u{0}ud+(x*), 

with B-(x*) = Xrjix*) nM^Q and B+(x*) = X^(x*) Since x* G B^(x*) and x* / 0, 

then it either belongs to B“(x*) or to D+(x;*). 

Suppose that x* G D“(x*), then we have that < 0 for every X = 1, - ■ ■ ,i. By Theorem 15.11 
every function Fa(xa;x*) ™ (l 66 ]l is monotonous decreasing. Thus, for every A, the following 
strict inequalities hold: 

Fa( 0 ;x*) <Fa(x^;X*). 

The summation over A results in: 

J;Fa(0;x*)< J;Fa(x1;x*) = 0, 

A A 

what is in contradiction with expression (j89p . A similar argument for x* G B+(x*) leads to: 

J;Fa(0;x*)>^Fa(x1;x*) = 0, 

A A 

which again is in contradiction with (|89li . This proves that for any fc G "H, the set of feasible 
solutions contains just one element f*, which corresponds to complex balanced equilibrium 
solutions. □ 

Proposition 6.2 For each k ^T-L, there exists exaetly one complex balanced equilibrium in each 
positive stoichiometrie compatibility class, and this equilibrium is locally asymptotically stable 
within the corresponding stoiehiometrie compatibility class. Furthermore, if the reaction network 
consists of one linkage class, then the asymptotic stability of the equilibrium point corresponding 
to any k within its stoichiometric compatibility class is global. 


Proof: As discussed above, for a given k £ T-L, the set of feasible solutions contains just one 
element f*. Let cq G be one (complex balanced) equilibrium point so that according to (IHSli 
we have that Inf* = S^lncQ. Then, any (complex balanced) equilibrium satisfies: 

S^(lnc — Inco) = 0, (91) 
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which coincides with the set IA{cq) defined as (jB.ip in Proposition IB. 11 From this proposition, 
it follows that the set contains exactly one element in each positive stoichiometric compatibility 
class. 

That each equilibrium is locally asymptotically stable follows from a standard result presented 
in [22] (Lecture 5, Proposition 5.3) and summarized in Appendix iBl fProposition IB.2p . The 
global stability of the complex balanced equilibria in the single linkage class case is proved in 

m- □ 

Remark 6.1 For weakly reversible deficiency zero reaction networks, Im(S''^) spans what 

in turn implies that any vector of reaction rate coefficients k G R>o will be an element of the 
Horn set Thus, from Propositions and \6.Sl it follows that any equilibrium solution 

will be unique (one in each positive stoichiometric compatibility class) and locally asymptotically 
stable. If the network is irreversible, positive equilibrium does not exist. Such conclusions have 
been formally stated in the so-called Defieiency Zero Theorem. I^ . Additionally, we remark 
here that aecording to recent results, the stability of any complex balanced equilibrium point is 
most probably global In- 

Remark 6.2 If the reaction network is elementary and compatible with thermodynamics (this 
implying detailed balancing), any allowed reaetion constant for the network must be in H. In 
this way, the definition of the Horn set could be considered as an alternative statement of the 
Wegscheider conditions (see for a classical statement of the conditions). 


Example: The Horn set for a weakly reversible reaction 


Let us consider the example discussed in Subsection 14.31 of a one linkage class network of defi¬ 
ciency <5 = 2. The corresponding graph structure, stoichiometry, and the set of possible param¬ 
eters is depicted in Figure [71 For this 2-parameter network, the Horn set (IHHIl is obtained by 
finding those (e,/3) that make InP orthogonal to and g^. Note that this implies that InP 
lies in Range(S'^). The conditions can be written as: 


1 0 -2 
0 1 -3 


In 


e//5 

1 

e//3 


= 0 . 


(92) 


This results in a set with parameters satisfying a = fi, that lead to complex balanced solu¬ 
tions. By construction, lnf(0;x) = 0 for every %, and therefore F(0;x) = 0. Since F(x;x) is 
monotonous decreasing, only complex balanced solutions exist for parameters in the Horn set. 


36 









6.2 The Deficiency One Theorem revisited 


Next, we present a result that might be a basis for an alternative proof of the deficiency one 
theorem for weakly reversible reaction networks. The theorem was originally proposed by [23] 
and recently discussed by [8], employing in both cases a graph theoretical formalism. The 
argument we propose builds on the following observations; On the one hand, the existence of a 
basis in the kernel of S that has at most one vector g'’ per linkage class. On the other hand, a 
particular orthogonal structure for such basis. 

The structure is such that for each element of the basis, the only possible nonzero coordinates 
must be at the location of the complexes that correspond to the linkage class the vector is 
associated to. Orthogonality on a normalized basis is formally expressed as: 

(g*)V=(5y, (93) 

where 6ij denotes the Kronecker delta. The structure of the vectors is illustrated in Figure 
El for a network consisting of 3 linkage classes, with the grey areas representing the non-zero 
vector coordinates. As sketched in the same figure, that particular structure of the g^ vectors 
decouples the corresponding feasibility conditions (Lemma 14.21) along linkage classes, so to have 
one feasibility function per linkage class. As we will see, monotonicity of such functions (Theorem 
EH) will ensure uniqueness. 


A=1 A=2 A=3 

I I I 


N 

—^ 

Fi(a;i) = 0 


g 


2 


g 


;? 


F2{x2) = 0 


F3(a^3) = 0 



Figure 9: Vectors g*' form a basis of the kernel of S where the grey areas indicate the only 
possible non-zero coordinates. This structure decouples feasibility conditions (at most one per 
linkage class). Monotonicity of Fa(xa), schematically represented as discontinuous lines at the 
right of the figure, then leads to just one solution per linkage class. 

Theorem 6.1 Let us consider a weakly reversible reaction network with i linkage classes such 
that 6 = where 6\ is either 0 or 1. Then, there will be a unique equilibrium in each 

positive stoichiometric compatibility class. 
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Proof: Let t be the number of linkage classes having deficiency 1. Since for each linkage class, 
5x can be either zero or 1, then t < i. Assume, without loss of generality, that the deficiencies of 
linkage classes A = 1,... ,t are 1, whereas for the remaining X = t + 1,... ,i, linkage deficiencies 
are <5^ = 0. 


For the first t linkage classes, let the non-zero vectors G IK'^A 1 (for A = 1,..., t) be the basis 
associated to the kernel of Sx so that 5 aP^ = 0. From the assumption, we have that: 


5 = = 

A 


and a basis {g^ | A = 1,..., t} for 
can be constructed as: 

isY = [ (pY 
{sY = [ ioY 
isY = [ {oY 


the kernel of S satisfying (|9c{p under proper normalization 


(Oa)^ ••• (0^)^ (0,f] 

(p^r ••• (0,)^ ••• (0,)^] (94) 

(Oi)^ ••• {pY ••• (O^)^] 


where Oa G for A = are zero vectors. Note that under the assumptions of the 

theorem, for any equivalent basis {g^ | A = 1,... ,t} of the kernel of S, the sub-vectors g^ in 
(flSl) must satisfy that ^AgA = 0 for A = 1,... , t. This is so since each element g'^ of the basis 
can be expressed as a linear combination of (j94p . what in addition implies that the vectors g^ 
and p''' are parallel. Equivalently, there exists a non-zero scalar //a such that g^ = ^aP"''- 


Using the basis {g^ | A = 1,... ,t}, conditions (|64p in Lemma 14.21 can be written as: 

(g^)^lnf^(x;x) = 0, for A = I,--- ,t. (95) 

Taking into account the structure of the basis components (I94p . conditions in (I95p reduce to: 

{pY In (^ + ziEfipi)= 0 

(p^)^ In (r^ + zaE-'p^) = 0 (96) 

(pT In {n + ztEi^p^)= 0 

where za = ^caXa for A = 1,... ,t. The left hand side of each of the expressions in ([96]) . that 
we denote as Fa(za) for (A = 1,... ,t), is of the form (17X1) (Section [5]). According to Theorem 
EH each of those functions is monotonous decreasing and because of (1721) . it becomes zero at a 
point G Xa. Therefore, the set of feasible solutions contains just one element: 

••• fxYx) ••• ifizt) (YY ••• inV], (97) 

with fA( 2 ;^) = fl + z^La for a = 1, • • • ,t. 
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Similarly to the proof of Proposition [621 let cq € M™q be one equilibrium point so that according 
to (|88l) . {zl,..., z*) = S’^ IncQ. Then, any equilibrium satisfies: 

5'^(ln c — In Co) = 0, 

which coincides with the set dEIl). Finally, from Proposition IB. 11 it follows that the set contains 
exactly one element in each positive stoichiometric compatibility class. □ 

As a final remark, we note that this result allows us to conclude uniqueness of equilibria in each 
positive stoichiometric compatibility class (although not stability), for networks with deficiency 
other than zero, provided that feasibility conditions in Lemma [4. 2 1 can be decoupled along linkage 
classes, as discussed above. 


6.3 A complex network satisfying the deficiency one theorem 


Let us consider a reaction network involving m = 7 chemical species, we label with capital letters 
from A to G. The (reversible) reaction steps that take place are: 


2A^B 
B^A + C 
C + E^2G 
A + D^ E 


BU2C 
2G^A + G 

E^E 


2C^D 

A + D^ E 


(98) 


This particular reaction network comprises n = 10 complexes and i = 3 linkage classes, we 
represent in graph form in Figure fTOl with explicit indication of the species iFigure fTOhl as well 
as in terms of numbered complexes (Figure [TUJd). 

(a) (b) 




Figure 10: Graph representation for the reaction network, (a) The species that are part of each 
complex are explicitly indicated, (b) The same graph described in terms of numbered complexes. 

Complexes are grouped by linkage class in the sets Ci = {1,2,3,4,5}, £2 = {6,7,8} and 
^3 = {Sj 10}. The stoichiometry associated to the complexes is given (column-wise) in the 
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following molecularity matrix fSection 12.ip : 


2 

0 

0 

0 

1 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

0 


( 99 ) 


Choosing ji = 1, j 2 = 6 and = 9, as the reference complexes, matrices Sx^ at the right of 
expression (IMD, become: 


■ -2 

-2 

-2 

-1 ■ 


■ -1 

-1 ■ 


■ 0 ■ 

1 

0 

0 

0 


0 

0 


0 

0 

2 

0 

1 


0 

0 


1 

0 

0 

1 

0 

, ^2 = 

-1 

-1 

, ^3 = 

0 

0 

0 

0 

0 


1 

0 


1 

0 

0 

0 

0 


0 

1 


0 

1- 

0 

0 

0 

0 . 


0 

0 . 


_ -2 _ 


Net reaction fluxes take the form: 

hifp) = + fc3,2V'3 + ^5,2'05 “ (^2,1 + ^2,3 + fc2,5)V’2 

(psW = k2,3-1^2 + fc4,3V'4 + h,3'lp5 “ {h,2 + ^3,4 + fc3,5)V’3 
(pii'tp) = - ki^slpi 

= k2,5'4’2 + k3^5ljj3 - (^ 5,2 + A:5,3)'!/’5 ( 100 ) 

(j)7{tp) = kQj^pe + ks,7^p8 - (fey,6 + fe7,8)V'7 
08(^) = fee,806 + fe7,807 - (fe8,6 + fe8,7)08 
0io(0) = fe9,lO09 - felO,901O 

The remaining fluxes 0i(0), 06(0) and 09(0), associated to the reference complexes, are ob¬ 
tained by means of relation ([8]). 


The dimension of the stoichiometric subspace, which coincides with the rank of matrix S = 
[ 5i S 2 53 ], is s = 6, and renders a network deficiency (5 = 10 — 3 — 6 = 1. Hence, the kernel 
of S is one dimensional with a basis g^ = (0 —1/2 0 1 0 0 0).As shown in Section 
12.21 we identify the following three sub-vectors in that solve (1161) : 

5I = ( 0 -1/2 0 1)^, 

= ( 0 0 )^, 

53 ' = (0)^ • 

The canonical representation of the equilibrium set will be expressed in terms of matrices M\ 
that appear in Eqn ()37p . From the expressions (llOOp for the fluxes, we have that: 


Ml 


-k 


1,2 

T,2 

0 

0 

0 


fe2,l 

— (fe2,l + fe2,3 + fe2,5) 
fe2,3 
0 

fe2,5 


0 0 0 

fe3,2 0 A;5^2 

— (fe3,2 + fe3,4 + fe3,5) fe4,3 fe5,3 

fe3,4 —fe4,3 0 


fe3,5 0 — (fe5,2 + fe5,3) 


( 101 ) 
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-(fee,7 + kQ^s) 

kjfi 

kg,6 


—kg, 10 kiQ,g 

kg, 10 —kio,g 

M 2 = 

ke,7 

— {krfi + kj^s) 

kg,7 

, M3 = 


^6,8 

^ 7,8 

-{kg,6 + kg,7) 



Comparing each matrix with the structure given in (1411) . we get for each linkage class: 


( 102 ) 


El = 


-{h,l + k2,3 + k2,5) 
k2,3 

0 

^2,5 


k3,2 

— (^3,2 + ks^i + /Cs^s) 
^3,4 
^3,5 


0 A:5,2 

^4,3 ^5,3 

-k4,3 0 

0 — (A:5,2 + ^ 5 , 3 ) 


(103) 


E 2 = 


— (^ 7 , 6 + ^ 7 , 3 ) ksj 

kjfi -{ksfi + ksj) 


T 

«1 = [ ^1,2 000] 02 = [ / c 6,7 ^6,8 

= [ ^ 2,1 000 ] 62 = [ ^ 7,6 ksfi 

Expressions of the form (|5ip . which describe the family of 


E 3 — [ —^ 10,9 ] 

] 03 = [ kg^io ] 

63 =[^ 10 , 9 ]^ 
solutions, become as follows: 


(104) 

(105) 


fi(xi) = +xihi 

f2(o:2) = ^ 

f3(x3) = ^3 

where vectors at the right hand side, for the parameters given in Table [H become: 

fj = (1.0000 30.1818 75.4545 2.1091)'^, = (0.2000 0.8571)^ 

hi = (-0.2500 - 7.5455 - 18.8636 - 1.0273)^ = (0.2500) 


(106) 


(107) 


As it can be seen in Figure fTTk. function Fi(xi) = (< 7 ^)^ Infi(xi) is monotonous decreasing, 


Table 1: Reaction rate coefficients for the network 


A:i ^2 = 2.0 A: 2 ,i = 2.0 ^ 2,3 = 16.0 /c 3^2 = 0.5 ^ 3^4 = 2.5 A: 4^3 = 1.0 

^2,5 = 1-2 A:5,2 = 1-0 /c3^5 = 0.1 k^^s = 1.0 kgj = 1.0 kj^Q = 1.0 

^7,8 = 10-0 ksj = 1.4 ksfi = 2.1 = 1.0 fcg^io = 1.0 kip^g = 4.0 


with one solution (intersection with the x-axis) at = —6.69 as asserted by Theorem 16.11 
Another example is represented in Figure fTTb . for a network with reaction rates as in Table [H 
except for constants A; 2 ,i and k ^^2 which now take the value 10. The domain of the function is 
now constrained to the interval Xi = with = —0.6154 and Lf = +0.3636. Vector 

hi for this parameter set becomes hi = (—0.0500 0.1857 0.1857 — 0.0786)^. For this case the 
function crosses the x-axis at x\ = —0.5854. 


In order to compute all possible equilibrium solutions in the concentration space, we make 
use of f^(x*) = [ ff (xi) ^ ], where ff (xi) for xj = —6.69 and xj = —0.5854 takes, 

respectively, the values: 


fi(xt) = ( 2.6725 80.6609 201.6523 8.9815 )^, 
fi(xi) = ( 0.2293 0.0056 0.0056 0.0746 f. 
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(a) 


(b) 




Figure 11: Function Fi(x) = ( 5 '^)^Infi(x) corresponding to linkage class 1. Note that for the 
other linkage classes the functions are zero since vectors and are identically zero. In (a) 
the function is represented for the set of rate coefficients in Table [TJ In (b) the function is 
represented for two modified rate constants, k 2 ,i = /c 5^2 = 10 . 0 . 

The set of equilibrium concentrations is then computed by solving Infj, = S'^ Inc. In particular, 
for this network we can express the first 6 chemical species in terms of species G, what leads 
to a straight line in the Inc-space, which intersects the interior of any positive stoichiometric 
compatibility class defined in (fT71) (with B = (^ 1/2 1 1/2 1 3/2 3/2 1 )^) in just one 

point. 


7 Conclusions 


This contribution concentrates on the study of feasibility conditions to identify admissible equi¬ 
libria for weakly reversible mass action law (MAT) systems. To that purpose, a flux-based form 
of the model equations describing the time evolution of the species concentration has been ex¬ 
ploited, in combination with results from the theory of linear compartmental systems to develop 
a canonical representation of the equilibrium set. Ingredients of such representation include the 
so-called family of solutions, with the corresponding positivity conditions, and the feasibility 
functions employed to characterize the set of feasible (equilibrium) solutions. 

One main result of this contribution is that the introduced feasibility functions are monotonously 
decreasing on their domain. This allows us to establish connections with classical results in 
CRNT related to the existence and uniqueness of equilibria within positive stoichiometric com¬ 
patibility classes. In particular, we employ monotonicity to identify regions in the set of possible 
reaction rate coefficients leading to complex balancing, and to conclude uniqueness of equilibria 
for a class of positive deficiency networks. It is our hope that the proposed results might sup¬ 
port the understanding of the deficiency one theorem from a different point of view, with the 
possibility of an alternative proof. 
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A number of examples of different complexity are employed to illustrate the notions presented 
and their relations. As the examples show, all components used for the characterization of 
equilibria, in particular the family of solutions and the feasibility functions, can be computed 
efficiently in an algorithmic way, even for large kinetic models. Future work will be focused on 
the constructive application of these functions for the computational search or design of networks 
with unique equilibria. 
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A Some results required to prove Theorem 15.1 


Lemma A.l Let H G be C-Metzler and such that: 


H — a, 


(A.l) 


with a > 0. Let p G M” (with p ^ 0) be a vector with m positive, r zero and n — m — r negative 
components satisfying: 


and 


Then: 


and 


Moreover: 


Vl> ■ ■ ■ >Pk> ■ ■ ■ >Pm> ^ > Pm+r+l > ' ' ' > Pi > ' ' ' > Pn, 
Pm+1 — ' ' ■ — Pm+r — 0 ; 


g = Hp. 


k 

''L2 di 0 for every 1 < k < m, 
i=l 


n 

> 0 for every m + r + l<£<n. 

i=l 


m 

gi < 0 and 

i=l 


9i>0. 

2=m+‘r+l 


(A.2) 

(A.3) 

(A.4) 

(A.5) 

(A.6) 


Proof: Multiplying both sides of (lA.ip by the scalar > 0 and subtracting the result from 
()A.3p . we get: 

H{p - pkln) = g + Pk^. (A.7) 

Summing the first k elements and reordering terms results in: 

k k-i / k \ n / k \ k 

X] (pj -Pk)+ Y1 X] (pj - Pk) + i-Pk) (^-8) 

i=l j=l Vi=l / j=k+l V*=l / *=1 

Since H is C-Metzler, according to Definition 13.11 for every j = 1,... ,k, with k = 1,...,re, we 
have that: 

k n 

^ij ~ ~^j ~ ^ij — 0 - 
i=l i=k+l 

The first term at the right hand side of (jA.Sp is non-positive since, by construction, pj — Pk ^ 0 
for j = 1,... ,k — l, and the above summations are non-positive. The second term is non-positive 
since for every j = k + 1,... ,n and i ^ j, Hij > 0 and pj —pk<0- Thus, relation (IA.4p follows, 
since a is a nonnegative vector and pk is positive for k = 1,... ,m, so the third term at the right 
hand side of (jA.811 is also non-positive. 
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In a similar way we prove (]A.5p . Substituting pi < 0 for in (jA.7p . we get: 


H{p - piln) = g + Pia. (A.9) 

Summing the elements of g from t' = m + r + l,...,n gives: 

^ ) (pj - Pi) + i-Pi) (A-10) 

i=e j=i+\ \i=l ) j=l \i=l ) i=£ 

Because H is C-Metzler, we have that: 

n 

Y any j = , n. 

i=i 

Thus, the first term at the right hand side of of ()A.10p is non-negative, since (j)j — pi) < 0 for 
j = £ + 1,..., n. The second term in the expression is also non-negative, since the off-diagonal 
elements of H are non-negative and {pj — pi) > 0 for j = 1, — 1. Finally, the last term in 

(|A.10I1 is non-negative due to the negativity of pi and the non-negativity of a. 

Strict inequalities (IA.6D can be proven in a straightforward manner from expressions (lA.Sp and 
(|A.10p . if the non-zero components of vector a > 0 are within the first m and last n — m — r 
entries. This would be the case, since the last terms at the right hand side in both equations 
would be strictly negative (with k = m), and positive (with ^ = m -|- r -|- 1), respectively. 


If the non-zero components are not within the first m, nor within the last n — m — r entries, the 
strict inequalities still hold. In order to prove this point, we express H as: 


H = 


Hu 

Hi2 ' 

H 21 

H 22 


(A.ll) 


where Hu G Let the first m components of vector a to be zero. Then, H 12 € 

in ()A.lip must necessarily have at least one positive element (any non-zero element must be 
positive because H is C-Metzler). Suppose, on the contrary, that H 12 is a zero matrix. Then, 
by using mh we have that: 


— 0 , 


which means that Hu, and consequently H, are not invertible, contradicting the fact that H 
is C-Metzler and therefore, non-singular. Since at least one entry of H 12 is positive, the second 
term at the right hand side of ()A.8h for k = m must be strictly negative. 


A similar line of arguments can be employed if the last n — m — r components of a are zero, with 
matrix H 22 G and H 21 G instead of Hu and H 12 . Now, we 

suppose that H 21 is a zero matrix, what combined with dATl leads to: 

H22^n—m—r — O' 
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Thus H 22 , and consequently H, are not invertible, what is in contradiction with the fact that 
H is C-Metzler and therefore, non-singular. Since at least one entry of H 21 must be positive, 
the second term at the right hand side of (jA.lOp . for i = n — m — r, must be strictly positive, 
completing the proof. □ 

Lemma A.2 Let X C K and consider the function G{x) : X 1 -^ R defined as: 

n 

G{x) = ^giQi{x), (A.12) 

i=l 

where gi are the coordinates of the vector g = Hp, with H and p as in Lemma \A.l[ For every 
i = 1,... ,n and x G X, let also have that: 

^ ^ Qkix) ^ ^ Qmix^ > 0 > Qm+r+li^') ^ ^ Ql{pf) ^ . 

Qm+l{x) = ■ ■ ■ = Qm+rix) = 0. 

Then, G{x) < 0 for every x G X. 

Proof: First, we note that (jA.12p can be re-written as: 

k m 

G{x) = {Qi{x)-Q 2 {x))gi+{Q 2 {x)-Q 2 ,{x)){gi+g 2 )^ - ^{Qk{x)-Qk+i{x))^gi^ - VQm{x)'^gi+ 

i=l i=l 

n n 

Qm+'r+l(^) ^ ^ 9i ^ ^ 9i ‘ ‘ ‘ (Qn(^) Qn—l(^))^n; 

2 =m+r+l i=£ 

where implicitly, each Qi is assumed to be a function of x. From ()A.13p . we have that Qfix) — 
Qj{x) > 0 for every Qi{x) > Qj{x) and x G X, what implies that {Qk{x) — Qk+i{x)) > 0 for 
every A: = 1,..., m — 1, and {Qfix) — Q(,-i{x)) < 0 for every ^ = m -|- r -|- 1,..., n. Thus, from 
Lemma lA.ll we have that: 

m n 

G (x) ^ ^ ^ 9j “t“ Qm-\-r-\-l (^) E 9j- 

j=l j=m+r+l 

The signs Qm{x) > 0, Qm+r+iix) < 0 as well as inequalities (|A.6p . from Lemma [A. 11 make the 
right hand side of the above expression strictly negative, what completes the proof. □ 

Proposition A.l Under the conditions of Lemma \A.1{ let the m positive and the n — m — r 
negative components of p £ R” satisfy: 

Pi > ■ ■ ■ > Ps-l > Ps > Ps+i > ■ ■ ■ > Pm > 0, (A.15) 

0 > Pm+r+l> ■■■ >Pt-l > Pt>Pt+l> ■■■ >Pn- (A.16) 

Then: 

s n 

'^gi<0 and ^^*>0, (A.17) 

2=1 i=t 

for some s = 1,...,m — 1 and t = m + r + 2,... ,n. 
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Proof: The line of arguments is similar to that employed in Lemma lA.ll to prove (IA.6I) . If 
the non-zero components of vector a > 0 are within the first s and the last n — t + 1 entries, 
it is straightforward to prove strict inequalities from expressions (lA.Sp and (|A.10p . for the last 
terms at the right hand side in both equations (with k = s and i = t) are strictly negative, and 
positive, respectively. 


If, on the other hand, the first s and the last n — t + 1 entries of a > 0 are zero, then matrix 
H, which can be expressed as in (lA.lip with Hu G must have for Hu G at least 

one positive element. Otherwise, from (jA.ip . we would have that: 

Hu^s = 0 , 

what contradicts the hypothesis that H is C-Metzler and therefore, non-singular. Since at least 
one entry of Hu must be positive, and because of (|A.15p . pj — ps > 0 (for j = s + 1,... ,n) the 
second term at the right hand side of ()A.8p . for k = s, is strictly negative, what proves the hrst 
inequality in (IA.17I) . 


In order to prove the second inequality, we make use of a similar argument with H 22 G 
in (lA.lip . to show that H 21 G must have at least one positive entry. Prom (lA.lbp . 

we also have that pj — pt > 0 (for j = 1,..., t — 1) so the second term at the right hand side of 
(lA.IOh . for £ = t, must be strictly positive. This proves the second inequality in (lA.lTh . □ 


Proposition A. 2 Let us consider the following set of ordered parameters pi > Pk ^ Pk+i > 
0 > pi-i > Pi > Pn, and functions Ilj : M —>■ M 0 / the form Lij{x) = ln(l -|- xpj), with pj being a 
given parameter within the ordered set. Then, we have that: 


lim ni(x) = lim n„(x) = —00, 

X+—>■—( — ) —( —) 

lim (ni(x) - nfc(x)) = lim (n„(x) - n£(x)) = -CX), 
lim (Ukix) - nfc+i(x)) = Cl, lim {Ui_i{x) - n^(x)) = C2, 

^P\' ^Pn ^ 

lim (nfc(x) - nfc+i(x)) = C3, lim (n£_i(x) - Iii{x)) = Ci, 

X—>-+oo 00 


(A.18) 

(A.19) 

(A.20) 

(A.21) 


where Ci, C 2 , C 3 and C 4 are constants. 


Proof: The limits in ()A.18p follow since Hi increases and decreases monotonically in their 
respective domains (—1/pi,-|-oo), (— 00 ,—1/pn)- In order to prove (|A.19p . we have that: 

lim nA:(x) = Uk{-l/pi), and lim n£(x) = Iii{-l/pn), 

^Pl' ^Pn ' 

which are (negative) constants, because pi > pk and pe > Pn ■ Using ()A.I8p . we then get (IA.19I1 . 
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In order to compute the limits in ()A.20p . we have that: 


0<1 — — <1 — Thus, lim (11^(2:) — nfc_|_i(x)) = In- fe/Pi) ^ 

Pl~ Pi ^ l-{Pk+l/pi)- 

Similarly: 

0 < 1 - ^ < 1 - Thus, lim (n^_i(x) - n^(x)) = In ^ ~ > q 

Pn~ Pn ’ ^ l-{Pi/Pn) ~ 

In proving (IA.21h . we have that: 

lim (nfc(a:)-nfc+i(x)) = In lim ^ ^ ^nd lim (n£_i(x)-n(>(x)) = In lim ^ 

a;—>-+oo X—>-+oo i -)- xp^j^i x—>■—oo x—>—oo i -|- xp£ 

Thus, by the theorem of I’Hopital, we have that: 

lim (nfc(g:) — nfc-i-i(x)) = In > 0, and lim (n£_i(x) — n£(a:)) = In < 0. 
x—>-+00 Pk +1 X—>■—oo p£ 


□ 


B Some convenient results on uniqueness and stability 


For the sake of completeness, here we summarize in the form of propositions, two fundamental 
results from CRNT on uniqueness and stability. The complete set of arguments can be found 
in [22]. 

Lemma B.l (see also [2j) LetV{x) : X —>■ M, with^ C its domain, a convex function with 
continuous derivatives in X, and v{x) : X ^ R"" he the gradient ofV{x). Then, the following 
inequalities hold for every x € X.- 


(i) {xi){x — xi) < V{x) — V{xi), for any xi € X. 

(ii) [n{x2) — n{xi)]^ {x2 — xi) > 0, for any xi, X 2 G X. 

Inequalities are strict whenever x ^ xi or xi ^ X2 in (i) and (ii), respectively. 

Proof: In order to prove the first part, choose any xi G X and construct a function Ri(x; xi) as 
the difference between V{x) and its supporting hyperplane at xi. The supporting hyperplane is 
of the form H{x;xi) = V (xi)+i^’^(xi)(x—xi), and Bi{x; xi) = P {x)—H{x; xi). By construction, 
the function is strictly positive, i.e. it is positive for all x G X other than xi, and result (i) follows 
in a straightforward manner, since i?i(x; xi) = V (x) — V (xi) — z/^(xi)(x — xi) > 0, what implies 
that V{x) — P(xi) > :^^(xi)(x — xi). 
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To prove the second part, we note that Bi{x\xi) is itself a convex function since S/xBi = 
v{x) — so its Hessian coincides with that of the convex function V{x). By using the same 

supporting hyperplane argument, we construct the following strictly positive definite function 
around some X 2 G X: 

B 2 {x; xi, X 2 ) = Bi{x; xi) - Hi(x 2 ; xi) - [z^(x 2 ) - (x - X 2 ) > 0, 

where the inequality holds for any x € X. In particular, it holds for x = xi, and therefore: 

Hi(x 2 ;xi) + [i^{x 2 ) - (xi - X 2 ) < 0, 

which implies that Hi(x 2 ;xi) < [i'{x 2 ) — i^(xi)]^ (x 2 — xi), and proves (ii). □ 

Proposition B.l (Corollary 4.14 f22]) Let cq € be a fixed reference. The set 

U{cq) = {c € M™q I 5"'"(lnc — In Co) = 0}, (B.l) 

contains exactly one element in each positive stoichiometric compatibility class. 

Proof: In proving uniqueness, suppose that there are two elements: c*, c** € IA{cq), that belong 
to the same stoichiometric compatibility class. Then, we have that ^^(Inc* — Inc**) = 0, what 
implies that (Inc* — Inc**) is orthogonal to the stoichiometric subspace H. Because c* and 
c** are assumed to be in the same compatibility class, the vector c* — c** must belong to the 
stoichiometric subspace H, and the following relation hold: 

(lnc*-lnc**)^(c*-c**) =0. (B.2) 

Using the convex function U(c) = c^(lnc — 1), with gradient i^(c) = Inc, and applying Lemma 
IB. II (condition (ii)), it follows that equality ()B.2p holds if and only if c* = c**, what proves that 
the set L/(co) can have at most one element in each positive stoichiometric compatibility class. 

As pointed out in [22], the question of existence (i.e. that each (positive) stoichiometric com¬ 
patibility class in fact meets W(co)) is somewhat more difficult to answer than uniqueness. The 
complete argument can be found in |22] (Proposition 4.13). □ 

Proposition B.2 (see also [22] ) Complex balanced equilibria are locally asymptotically stable 
in all positive stoichiometric compatibility classes. 


Proof: First of all, let us make use of Eqn ([5|) to write the right hand side of system (j3|) as a 
summation over A, of functions: 

■ iVj - Vi)- (B-3) 

i&Cx j&Xi 
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Select some positive reference c* > 0 (its associated vector -0* is strictly positive) and re-write 
the previous expression in the equivalent form: 

• iVj - Vi)^ (B-4) 

ieCx j&Ti 

where 17 = Inc — Inc*. The inner product between V and TZ^{c) (IB.4p results into the following 
scalar function: 

V^n^{V) = Y “ ^*(^))’ (^' 5 ) 

i&Cx jeli 

where Zi{Ty) = yfv. In order to get an upper bound for (|B.5p . we make use of Lemma IB. II 
(condition (i)), with the convex function V{z) = e^, to obtain: 


e^‘(z,-z,)<e^^-e^^ 


(B.6) 


For any scalars Zi and Zj. Strict convexity of V(z) ensures that the equality holds only if Zi = zj. 
We also have that: 

n 

e-. -e^* = (£,-£,)^^£fce^F (B.7) 

k=l 

Combining (|B.7I) with (IB.6I1 . and substituting the resulting expression in (|B.5p . we get: 


v^TZ^{v) < 




Aim- 


Y 

i€Cx j&Ii 



(B.8) 


If the reference corresponds with a complex balanced equilibrium, then for every X = 1,... ,i, 
= 0 s-iid so is the right hand side of (IB.81) . Note that inequality is strict, in the sense 
that it holds whenever Zi 0 Zj, for every i,j G Cx. 


Local asymptotic stability is proved by the standard Lyapunov stability method (see for instance 
[39]) with the following Lyapunov function candidate, constructed as in the proof of Lemma fB.il 

B(c; c*) = V{c) - V{c*) - i/'^(c*)(c - c*) > 0, 

with 14(c), being a convex function of the form: 

m 

14(c) = ^Cj(lncj - 1). 

i=l 

Computing the derivative of B along dH), and using (jB.Sp . we get: 

i / n \ ^ 

B = ^I7^7^"(I7) < Y^^m = 0- (B.9) 

A=1 \i=l / A=1 

The result then follows, since B(c; c*) > 0 and B(c; c*) < 0, with equality only if c = c*. □ 
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